#######################################
### LIEBERMAN, MARTIN, MCMURRY 2021 ###
######### FIGURES AND TABLES ##########
#######################################

## Packages
library(foreign)
library(readstata13)
library(MASS)
library(multiwayvcov)
library(sandwich)
library(lmtest)
library(survival)
library(boot)
library(pastecs)
library(xtable)
library(stargazer)
library(plyr)
library(AER)
library(sensemakr)
# for mapping
library(rgdal)
library(plyr)
library(ggplot2)
library(rgeos)
library(maptools)
library(sp)
library(sf)
library(RColorBrewer)
library(broom)
library(ggmap)
library(tmap)
library(grid)
library(gridBase)
library(xtable)
library(sensemakr)

### LOAD DATA ###

# Gauteng Ward-level
gcro <- read.csv("gautengMerged.csv")

# Subset to wards without by-elections
gcro_noBE <- gcro[which(gcro$BE == 0),]
# Create subsets for analysis
gcro_noMV <- gcro_noBE[-which(gcro$Municipality == "GT422 - Midvaal [Meyerton]"), ] # excluding Midvaal
gcro_ANC <- gcro_noBE[which(gcro_noBE$Party == "AFRICAN NATIONAL CONGRESS"), ]
gcro_ANC_muns <- gcro_ANC[-which(gcro_ANC$Municipality == "GT422 - Midvaal [Meyerton]"), ] # excluding Midvaal, only non-ANC municipality in GT
gcro_DA <- gcro_noBE[which(gcro$Party == "DEMOCRATIC ALLIANCE"), ]
gcro_DA_noMV <- gcro_DA[-which(gcro_DA$Municipality == "GT422 - Midvaal [Meyerton]"), ] # excluding Midvaal
# Non-Fully Serviced Wards (less than 90% for any service)
gcro_nfs <- gcro_noBE[-which(gcro_noBE$water_source_13 > 0.9 & gcro_noBE$elec_13 > 0.9 & gcro_noBE$flush_toilet_13 > 0.9 & gcro_noBE$waste_disposal_13 > 0.9), ]
# Subset this to ANC only - not enough observations to subset to DA
gcro_nfs_ANC <- gcro_nfs[which(gcro_nfs$Party.Name == "AFRICAN NATIONAL CONGRESS"), ]

# National Municipal-level
d <- read.csv("munMerged.csv")
# Create subsets for analysis
d_ward <- d[d$Ward == 1, ]
d_pr <- d[d$PR == 1, ]
drule <- d[which(d$rulingparty_member == 1), ]
drule_ward <- drule[which(drule$Ward == 1), ]
drule_pr <- drule[which(drule$PR == 1), ]
drule_ANC <- drule[which(drule$winner == "AFRICAN NATIONAL CONGRESS"), ]
dnonrule <- d[which(d$rulingparty_member == 0), ]
# Only top 10 PR councillors
drule_top10 <- drule[drule$pr_listorder %in% c(1:10), ]
drule_ANC_top10 <- drule_ANC[drule_ANC$pr_listorder %in% c(1:10), ]
dnonrule_top10 <- dnonrule[dnonrule$pr_listorder %in% c(1:10), ]

# Afrobarometer Individual-Level (R6)
ab <- read.csv("afrobarometer_r6.csv")

# Shapefiles (for maps)
# provinces
provs <- readOGR(dsn = "./PR_SA_2011", layer = "PR_SA_2011")
# munis
munis <- readOGR(dsn = "./LocalMunicipalities2011", layer = "LocalMunicipalities2011")
# wards
wards <- readOGR(dsn = "./Wards2011", layer = "Wards2011")
# subset to gauteng
gt_wards <- wards[wards@data$PROVINCE == "Gauteng", ]
rm(wards)
# municipal-level covarates for map and appendix figure 
d_muni <- read.csv("muncovars.csv")

#### TABLES ####

## Table 1: Determinants of Councillor Renomination in Gauteng (Logit) ##

# Column 1: All Wards - satisfaction only
all.mod <- glm(renom ~ 
                 local_satisfied_15 
               + high_income_avg 
               + post_secondary_avg
               + black_avg 
               + ce_index15 
               + logpop2011 
               + win_margin2011
               + years_incumbent 
               + ever_switchedparty 
               + I(Municipality),
               data = gcro_noBE, family = binomial(link = "logit"))
# Column 2: All Wards - councillor satisfaction and service satisfaction
all.mod.subj <- glm(renom ~ 
                      local_satisfied_15 
                    +  subj_sp_index15
                    + high_income_avg 
                    + post_secondary_avg
                    + black_avg 
                    + ce_index15 
                    + logpop2011 
                    + win_margin2011
                    + years_incumbent 
                    + ever_switchedparty
                    + I(Municipality),
                    data = gcro_noBE, family = binomial(link = "logit"))
# Column 3: All wards - councillor satisfaction, service satisfaction, objective services
all.mod.obj <- glm(renom ~ 
                     local_satisfied_15 
                   + subj_sp_index15
                   + obj_sp_delta_index_1315
                   + formal_housing_delta_1315 
                   + unemp_delta_1315_reverse
                   + high_income_avg 
                   + post_secondary_avg
                   + black_avg 
                   + ce_index15 
                   + logpop2011 
                   + win_margin2011
                   + years_incumbent 
                   + ever_switchedparty 
                   + I(Municipality),
                   data = gcro_noBE, family = binomial(link = "logit"))
# Column 4: Councillor satisfaction only, ANC wards only
anc.mod <- glm(renom ~ 
                 local_satisfied_15
               + high_income_avg
               + post_secondary_avg
               + black_avg 
               + ce_index15 
               + logpop2011 
               + win_margin2011
               + years_incumbent 
               + ever_switchedparty
               + I(Municipality),
               data = gcro_ANC, family = binomial(link = "logit"))
# Column 5: councillor satisfaction and service satisfaction, ANC only
anc.mod.subj <- glm(renom ~ 
                      local_satisfied_15 
                    + subj_sp_index15
                    + high_income_avg 
                    + post_secondary_avg
                    + black_avg 
                    + ce_index15 
                    + logpop2011 
                    + win_margin2011
                    + years_incumbent 
                    + ever_switchedparty 
                    + I(Municipality),
                    data = gcro_ANC, family = binomial(link = "logit"))
# Column 6: councillor satisfaction, service satisfaction, and objective sp, ANC only
anc.mod.obj <- glm(renom ~ 
                     local_satisfied_15 
                   + subj_sp_index15 
                   + obj_sp_delta_index_1315 
                   + formal_housing_delta_1315 
                   + unemp_delta_1315_reverse
                   + high_income_avg 
                   + post_secondary_avg
                   + black_avg
                   + ce_index15
                   + logpop2011
                   + win_margin2011
                   + years_incumbent
                   + ever_switchedparty 
                   + I(Municipality),
                   data = gcro_ANC, family = binomial(link = "logit"))
# Column 7: councillor satisfaction only, DA wards
da.mod <- glm(renom ~ 
                local_satisfied_15
              + high_income_avg
              + post_secondary_avg
              + black_avg
              + ce_index15
              + logpop2011
              + win_margin2011
              + years_incumbent
              + ever_switchedparty 
              + I(Municipality),
              data = gcro_DA, family = binomial(link = "logit"))
# Column 8: councillor satisfaction and service satisfaction, DA only
da.mod.subj <- glm(renom ~ 
                     local_satisfied_15
                   + subj_sp_index15
                   + high_income_avg
                   + post_secondary_avg
                   + black_avg
                   + ce_index15
                   + logpop2011
                   + win_margin2011
                   + years_incumbent
                   + ever_switchedparty 
                   + I(Municipality),
                   data = gcro_DA, family = binomial(link = "logit"))
# Column 9: councillor satisfaction, service satisfaction, and objective SP, DA only
da.mod.obj <- glm(renom ~ 
                    local_satisfied_15
                  + subj_sp_index15
                  + obj_sp_delta_index_1315
                  + formal_housing_delta_1315
                  + unemp_delta_1315_reverse
                  + high_income_avg
                  + post_secondary_avg
                  + black_avg
                  + ce_index15
                  + logpop2011
                  + win_margin2011
                  + years_incumbent
                  + ever_switchedparty 
                  + I(Municipality),
                  data = gcro_DA, family = binomial(link = "logit"))
# Create table
labs <- c("Councillor Satisfaction", "Service Satisfaction Index",
          "Service Coverage $\\Delta$", "Formal Housing $\\Delta$",
          "Employment $\\Delta$", "High Income (\\%)", "Post Secondary (\\%)",
          "Ethnicity - African (\\%)", "Civic Engagement",
          "Log Population (2011)", "Win Margin (2011)", "Years Incumbent", "Switched Party")
stargazer(all.mod,
          all.mod.subj,
          all.mod.obj,
          anc.mod,
          anc.mod.subj,
          anc.mod.obj,
          da.mod,
          da.mod.subj,
          da.mod.obj,
          omit = c("I\\(", "Constant"),
          no.space = TRUE,
          covariate.labels = labs,
          column.labels = c(rep("All", 3), rep("ANC", 3), rep("DA", 3)),
          omit.stat = c("ll", "aic"),
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = F,
          font.size = "footnotesize",
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes = "*$p<0.05$",
          notes.label = "",
          title = "Determinants of Councillor Renomination in Gauteng (Logit)",
          add.lines = list(c("Municipality FE", rep("Yes", 9))),
          out = "wardmain.tex",
          label = "tab:wardmain"
)

## Ward Main Effects simulation - substantive effect size in text ##
mycoefs <- coef(all.mod.obj)
set.seed(02139)
simbetas <- mvrnorm(10000, mycoefs, vcov(all.mod.obj))
medians <- apply(model.matrix(all.mod.obj), 2, median)
myvars <- c("local_satisfied_15", "subj_sp_index15",
            "obj_sp_delta_index_1315",
            "formal_housing_delta_1315", "unemp_delta_1315_reverse")
allmod.res <- matrix(nrow = 3, ncol = length(myvars))
for(i in 1:length(myvars)){
  var.sd <- sd(na.omit(gcro_noBE[, myvars[i]]))
  medians_hi <- medians
  medians_hi[myvars[i]] <- medians[myvars[i]] + var.sd
  pps_med <- inv.logit(simbetas %*% medians)
  pps_hi <- inv.logit(simbetas %*% medians_hi)
  pps_diff <- pps_hi - pps_med
  allmod.res[, i] <- quantile(pps_diff, c(0.025, 0.5, 0.975))
}
allmod.res

## Table 2: Determinants of Councillor Renomination - All Municipalities (Logit) ## 

# Column 1: All ruling party councillors, subjective ratings only
munmod.all <- glm(renom ~ 
                  + service_rating_index2
                  + logpopchange
                  + black_avg
                  + post_secondary_avg
                  + win_margin
                  + Ward
                  + years_incumbent
                  + ever_switchedparty 
                  + I(Province)
                  ,
                  data = drule, family = binomial(link = "logit"))
munmod.all.c <- coeftest(munmod.all, vcov = cluster.vcov(munmod.all, drule$cat_b))
# Column 2: ANC only, subjective ratings only
munmod.alla <- glm(renom ~
                     service_rating_index2
                   + logpopchange
                   + black_avg
                   + post_secondary_avg
                   + win_margin
                   + Ward
                   + years_incumbent
                   + ever_switchedparty 
                   + I(Province)
             ,
                   data = drule_ANC, family = binomial(link = "logit"))
munmod.alla.c <- coeftest(munmod.alla, vcov = cluster.vcov(munmod.alla, drule_ANC$cat_b))
# Column 3: Non-ruling party councillors only, subjective ratings only
munmod.nonrule <- glm(renom ~
                        service_rating_index2
                      + logpopchange
                      + black_avg
                      + post_secondary_avg
                      + win_margin
                      + Ward
                      + years_incumbent
                      + ever_switchedparty 
                      + I(Province)
                 , 
                      data = dnonrule, family = binomial(link = "logit"))
munmod.nonrule.c <- coeftest(munmod.nonrule, vcov = cluster.vcov(munmod.nonrule, dnonrule$cat_b))
# Column 4: All ruling party councillors, controlling for objective service improvements
munmod.all.obj <- glm(renom ~
                        service_rating_index2
                      + service_pct_change_index1
                      + formaldwell_pct_delta
                      + logpopchange
                      + black_avg
                      + post_secondary_avg
                      + win_margin
                      + Ward
                      + years_incumbent
                      + ever_switchedparty
                      + I(Province),
                      data = drule, family = binomial(link = "logit"))
munmod.all.obj.c <- coeftest(munmod.all.obj, vcov = cluster.vcov(munmod.all.obj, drule$cat_b))
# Column 5: ANC councillors only, controlling for objective service improvements
munmod.alla.obj <- glm(renom ~  
                      service_rating_index2
                      + service_pct_change_index1
                      + formaldwell_pct_delta
                      + logpopchange
                      + black_avg
                      + post_secondary_avg
                      + win_margin
                      + Ward
                      + years_incumbent
                      + ever_switchedparty
                      + I(Province),
                      data = drule_ANC, family = binomial(link = "logit"))
munmod.alla.obj.c <- coeftest(munmod.alla.obj, vcov = cluster.vcov(munmod.alla.obj, drule_ANC$cat_b))
# Column 6: Non-ruling party councillors only, controlling for objective service improvements
munmod.nonrule.obj <- glm(renom ~
                            service_rating_index2
                          + service_pct_change_index1
                          + formaldwell_pct_delta
                          + logpopchange
                          + black_avg
                          + post_secondary_avg
                          + win_margin
                          + Ward
                          + years_incumbent
                          + ever_switchedparty
                          + I(Province),
                          data = dnonrule, family = binomial(link = "logit"))
munmod.nonrule.obj.c <- coeftest(munmod.nonrule.obj, vcov = cluster.vcov(munmod.nonrule.obj, dnonrule$cat_b))
# Create table
labs <- c("Service Rating Index", "Service Coverage $\\Delta$", 
          "Formal Housing $\\Delta$", "Log Population $\\Delta$",
          "Ethnicity - African (\\%)", "Post Secondary (\\%)", "Win Margin (2011)", "
          Ward Councillor",
          "Years Incumbent", "Switched Party")
nmunis <- length(unique(drule$cat_b))
nmunis_anc <- length(unique(drule_ANC$cat_b))
stargazer(munmod.all.c,
          munmod.alla.c,
          munmod.nonrule.c,
          munmod.all.obj.c,
          munmod.alla.obj.c,
          munmod.nonrule.obj.c,
          no.space = TRUE,
          covariate.labels = labs,
          column.labels = rep(c("All Ruling", "ANC", "Non-Ruling"), 2),
          omit = c("Province", "Constant"),
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = FALSE,
          add.lines = list(
            c("N Councillors", nobs(munmod.all), nobs(munmod.alla), nobs(munmod.nonrule), nobs(munmod.all.obj),
              nobs(munmod.alla.obj), nobs(munmod.nonrule.obj)),
            c("N Municipalities", rep(c(nmunis, nmunis_anc, nmunis), 2)),
            c("Province FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")
          ),
          font.size = "footnotesize",
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes = "*$p < 0.05$",
          notes.label = "RSE clustered at municipal level",
          title = "Determinants of Councillor Renomination - All Municipalities (Logit)",
          label = "tab:munmainupdated",
          out = "munmain_updated.tex"
         )

## MUNI MAIN EFFECTS SIMULATION- FOR NUMBERS IN TEXT (not tables) ##
mainmunFun <- function(mod, data){
  mycoefs <- coef(mod)
  set.seed(02139)
  simbetas <- mvrnorm(10000, mycoefs, cluster.vcov(mod, as.character(data$cat_b)))
  medians <- apply(model.matrix(mod), 2, median)
  vars <- c("service_rating_index2", "service_pct_change_index1", "formaldwell_pct_delta")
  res <- matrix(nrow = 3, ncol = length(vars))
  for(i in 1:length(vars)){
    var.sd <- sd(na.omit(data[, vars[i]]))
    medians_hi <- medians
    medians_hi[vars[i]] <- medians[vars[i]] + var.sd
    pps_med <- inv.logit(simbetas %*% medians)
    pps_hi <- inv.logit(simbetas %*% medians_hi)
    pps_diff <- pps_hi - pps_med
    res[, i] <- quantile(pps_diff, c(0.025, 0.5, 0.975))
  }
  return(res)
}
munres_allrule <- mainmunFun(munmod.all.obj, drule)
munres_allrule

### FIGURES ###

## Figure 1: Distribution of nomination outcome (Muni) ##

table(d$renom_type) / nrow(d)

pdf("promotion_hist_muni_comb.pdf")
par(mfrow = c(3, 1))
xx <- barplot(table(d$renom_type_num),
              names.arg = c("no nomination", "ward \n nomination",
                            "PR \n nomination", "nominated to \n higher office"),
              ylim = c(0, 4500),
              main = "2016 Nomination Outcomes \n All 2011 Councillors (n = 8,377)",
              cex.main = 0.9,
              cex.names = 0.7)
text(x = xx,
     y = as.numeric(table(d$renom_type_num)),
     labels = as.numeric(table(d$renom_type_num)),
     pos = 3,
     cex = 0.8)
xx <- barplot(table(d_ward$renom_type_num),
              names.arg = c("no nomination", "ward \n nomination",
                            "PR \n nomination", "nominated to \n higher office"),
              ylim = c(0, 4500),
              #main = "2016 Nomination Outcomes \n 2011 Ruling Party Councillors (All)",
              main = "2011 Ward Councillors (n = 4,220)",
              cex.main = 0.9,
              cex.names = 0.7)
text(x = xx,
     y = as.numeric(table(d_ward$renom_type_num)),
     labels = as.numeric(table(d_ward$renom_type_num)),
     pos = 3,
     cex = 0.8)
xx <- barplot(table(d_pr$renom_type_num),
              names.arg = c("no nomination", "ward \n nomination",
                            "PR \n nomination", "nominated to \n higher office"),
              ylim = c(0, 4500),
              main = "2011 PR Councillors (n = 4,157)",
              cex.main = 0.9,
              cex.names = 0.7)
text(x = xx,
     y = as.numeric(table(d_pr$renom_type_num)),
     labels = as.numeric(table(d_pr$renom_type_num)),
     pos = 3,
     cex = 0.8)
dev.off()

## Figure 2: Satisfaction with local ward councillor in Gauteng wards ##

# Merge in covariates
gt_wards@data$WARD_ID_NUM <- as.numeric(as.character(gt_wards@data$WARD_ID))
gt_wards <- merge(gt_wards, gcro, by.x = "WARD_ID_NUM", by.y = "ward", all.x = TRUE)
# Create color palette
pal_greys <- brewer.pal(7, "Greys")

pdf("map_ward_lcsatisfaction_paper_inlay.pdf")
plot.new()
vp_1 <- viewport(x = 0,
                 y = 0,
                 width = 0.91,
                 height = 1,
                 just = c("left", "bottom"))
vp_2 <- viewport(x = 0.21,
                 y = 0.8,
                 width = 0.22,
                 height = 0.25,
                 just = c("centre"))
vp_3 <- viewport(x = 0.21,
                 y = 0.8,
                 width = 0.23,
                 height = 0.26,
                 just = c("centre"))
# main map
pushViewport(vp_1)
par(new = TRUE, fig = gridFIG())
spplot(gt_wards, "local_satisfied_15", col.regions = pal_greys, cuts = 5,
       par.settings = list(axis.line = list(col = 'transparent'))
       #,
       #main = "Average Local Councillor Satisfaction Rating 2015 \n Gauteng Province"
)
# inset map
pushViewport(vp_3)
par(new = TRUE, fig = gridFIG(), mar = rep(0, 4))
grid.rect()
grid.rect(width=unit(1, "npc")-unit(0.03, 'lines'),
          height = unit(1, "npc") - unit(0.03, 'lines'))

upViewport()
pushViewport(vp_2)
par(new = TRUE, fig = gridFIG(), mar = rep(0, 4))
plot(provs)
plot(provs[provs$PR_MDB_C == "GT", ], col = pal_greys[7], add = T)
dev.off()

## Figure 3: Index of satisfaction with public services in South African municipalities ##

# Merge in covariates
munis <- merge(munis, d_muni, by.x = "CAT_B", by.y = "cat_b", all.x = TRUE)
# Create color palette
pal_greys <- brewer.pal(7, "Greys")

pdf("map_mun_servicerating_paper.pdf")
spplot(munis, "service_rating_index2", col.regions = pal_greys, cuts = 6,
       #main = "Service Rating Index 2015 (All Municipalities)",
       par.settings = list(axis.line = list(col = 'transparent')))
dev.off()

## Figure 4: Effects of Councillor Satisfaction on Predicted Probability of Councillor Renomination (Ward-Level) ##

lc.modx <- glm(renom ~
                 local_satisfied_15*win_margin2011
               + obj_sp_delta_index_1315
               + obj_sp_index13
               + formal_housing_delta_1315
               + unemp_delta_1315_reverse
               + subj_sp_index15 
               + high_income_avg
               + post_secondary_avg
               + black_avg
               + ce_index15
               + logpop2011
               + years_incumbent
               + ever_switchedparty 
               + I(Municipality),
               data = gcro_noBE, family = binomial(link = "logit"))
mycoefs <- coef(lc.modx)
simbetas <- mvrnorm(10000, mycoefs, vcov(lc.modx))
medians <- apply(model.matrix(lc.modx), 2, median)
simFunX <- function(mod, sd.var, winmargin){
  mycoefs <- coef(mod)
  set.seed(02139)
  simbetas <- mvrnorm(10000, mycoefs, vcov(mod))
  medians <- apply(model.matrix(mod), 2, median)
  medians[3] <- winmargin
  numcoefs <- length(medians)
  medians[numcoefs] <- medians[2] * winmargin
  medians_hi <- medians
  medians_hi[2] <- medians[2] + sd.var
  medians_hi[numcoefs] <- medians_hi[2] * winmargin
  pps_med <- inv.logit(simbetas %*% medians)
  pps_hi <- inv.logit(simbetas %*% medians_hi)
  fdiffs <- pps_hi - pps_med
  res <- quantile(fdiffs, c(0.025, 0.5, 0.975))
  return(res)
}
winmargins <- seq(0, 1, 0.1)
simRes <- sapply(winmargins, simFunX, mod = lc.modx, sd.var = sd(na.omit(gcro_noBE$local_satisfied_15)))
simRes

par(mar = c(5.1, 4.1, 4.1, 2.1))
pdf("lc_satisfaction_competition_notitle.pdf")
plot(x = winmargins,
     y = simRes[2, ],
     ylim = c(-0.15, 0.35),
     pch = 20,
     xlab = "Win Margin (2011)",
     ylab = "Change in Predicted Prob. of Renomination",
     #main = "Local Councillor Satisfaction and Renomination Probability \n Gauteng Wards",
     cex.axis = 0.7,
     cex.main = 0.9)
abline(h = 0, lty = 3)
segments(x0 = winmargins,
         x1 = winmargins,
         y0 = simRes[1, ],
         y1 = simRes[3, ])
rug(gcro_noBE$win_margin2011)
dev.off()

## Figure 5: Service satisfaction and PR promotion outcomes (Top 10 only)

ord.mod2 <- polr(factor(newprom) ~ service_rating_index2 + service_pct_change_index1 + formaldwell_pct_delta + logpopchange + black_avg + post_secondary_avg + win_margin + years_incumbent + ever_switchedparty + position_2011, data = drule_top10, method = "logistic")
ord.mod2.FE <- polr(factor(newprom) ~ service_rating_index2 + service_pct_change_index1 + formaldwell_pct_delta + logpopchange + black_avg + post_secondary_avg + win_margin + years_incumbent + ever_switchedparty + position_2011 + I(Province), data = drule_top10, method = "logistic")
betas <- ord.mod2.FE$coefficients
var_betas <- vcov(ord.mod2.FE)

set.seed(02139)
simbetas <- mvrnorm(10000, betas, var_betas[1:length(betas),1:length(betas)])
y <- as.matrix(drule_top10$newprom)
# Service Ratings
service_value1 <- median(drule_top10$service_rating_index2) 
service_value2 <- median(drule_top10$service_rating_index2) + sd(drule_top10$service_rating_index2)
y.0.effect <- c()
y.1.effect <- c()
y.2.effect <- c()
for(j in 1:nrow(simbetas)){
  mm <- model.matrix(ord.mod2.FE)
  mm <- apply(mm,2,median)
  mm[2] <- service_value1
  XB = simbetas[j, ]*mm[2:(length(betas)+1)]
  psi.0 <- 0
  psi.1 <- 1 #create psi parameters
  psi.2 <- 2
  #calculate predicted probabilities
  result <- c()
  for(i in 1:length(y)){
    #if y = 0
    if(y[i]==0){result[i] <- sum(log((exp(psi.0 - XB)) / (1 + exp(psi.0 - XB))))}
    #if y = 1
    else if(y[i]==1){result[i] <- sum(log(((exp(psi.1 - XB)) / (1 + exp(psi.1 - XB))) - ((exp(psi.0 - XB)) / (1 + exp(psi.0 - XB)))))}
    #if y = 2
    else if(y[i]==2){result[i] <- sum(log(((exp(psi.2 - XB)) / (1 + exp(psi.2 - XB))) - ((exp(psi.1 - XB)) / (1 + exp(psi.1 - XB)))))}
  }
  mm2 <- model.matrix(ord.mod2.FE)
  mm2 <- apply(mm2, 2, median)
  mm2[2] <- service_value2
  XB2 <- simbetas[j,]*mm2[2:(length(betas)+1)]
  result2 <- c()
  for(i in 1:length(y)){
    #if y = 0
    if(y[i]==0){result2[i] <- sum(log((exp(psi.0 - XB2)) / (1 + exp(psi.0 - XB2))))}
    #if y = 1
    else if(y[i]==1){result2[i] <- sum(log(((exp(psi.1 - XB2)) / (1 + exp(psi.1 - XB2))) - ((exp(psi.0 - XB2)) / (1 + exp(psi.0 - XB2)))))}
    #if y = 2
    else if(y[i]==2){result2[i] <- sum(log(((exp(psi.2 - XB2)) / (1 + exp(psi.2 - XB2))) - ((exp(psi.1 - XB2)) / (1 + exp(psi.1 - XB2)))))}
  }
  diff = result2 - result
  #average difference over categories
  y.0.effect[j] <- mean(diff[which(y==0)])
  y.1.effect[j] <- mean(diff[which(y==1)])
  y.2.effect[j] <- mean(diff[which(y==2)])
}
# Demotion - service rating
mean(y.0.effect)
quantile(y.0.effect, c(0.025, 0.975))
quantile(y.0.effect, c(0.05, 0.95))

# Same position - service rating
mean(y.1.effect)
quantile(y.1.effect, c(0.025, 0.975))
quantile(y.1.effect, c(0.05, 0.95))

# Promotion - service rating
mean(y.2.effect)
quantile(y.2.effect, c(0.025, 0.975))
quantile(y.2.effect, c(0.05, 0.95))

pdf("orderedlogit_service_rating_top10.pdf")
par(mfrow = c(1,1))
plot(y = mean(y.0.effect), x = 1, ylim = c(-.2, .2), xlim = c(.9, 1.5), pch = 20,
     xaxt = 'n', ylab = "Change in Pr Probability (Y|X)", xlab = " ",
     # main = "Service Rating Increase", cex.lab=1.2)
     cex.lab=1.2)
quantiles0 = quantile(y.0.effect, c(0.025, 0.975))
segments(x0 = 1, x1 = 1, y0 = quantiles0[1], y1 = quantiles0[2])
points(y = mean(y.1.effect), x = 1.2, pch = 20)
quantiles1 = quantile(y.1.effect, c(0.025, 0.975))
segments(x0 = 1.2, x1 = 1.2, y0 = quantiles1[1], y1 = quantiles1[2])
points(y = mean(y.2.effect), x = 1.4, pch = 20)
quantiles2 = quantile(y.2.effect, c(0.025, 0.975))
segments(x0 = 1.4, x1 = 1.4, y0 = quantiles2[1], y1 = quantiles2[2])
abline(h = 0)
par(xpd = NA)
text(x = 1, y = -0.24, "Demoted", cex = 1)
text(x = 1.2, y = -0.24, "Same Position/Ward", cex = 1)
text(x = 1.4, y = -0.24, "Promoted", cex = 1)
dev.off()


#### APPENDIX TABLES ####

## Appendix Table 1: Summary Statistics for Ward-Level Analysis (Gauteng Province) ##

gsum <- gcro_noBE[, c("renom",
                      "years_incumbent", "ever_switchedparty",
                      "obj_sp_delta_index_1315", "local_satisfied_15",
                      "subj_sp_index15", "formal_housing_delta_1315",
                      "unemp_delta_1315_reverse", "high_income_avg",
                      "post_secondary_avg", "black_avg", "ce_index15",
                      "logpop2011",
                      "win_margin2011")]
sumtab_ward <- t(stat.desc(gsum))[, c("nbr.val", "min", "max", "range", "median",
                                      "mean", "std.dev")]
varnames <- c("Renomination",
              "Years Incumbent", "Switched Party", 
              "Service Coverage $\\Delta$ (2013-2015)",
              "Councillor Satisfaction", "Service Satisfaction Index",
              "Formal Housing $\\Delta$ (2013-2015)", "Employment $\\Delta$ (2013-2015)",
              "High Income (\\%)", "Post Secondary (\\%)", "Ethnicity - African (\\%)",
              "Civic Engagement", "Log Population (2011)",
              "Win Margin (2011)")
row.names(sumtab_ward) <- varnames
colnames(sumtab_ward) <- c("N obs", "Min", "Max", "Range", "Median", "Mean", "Std. Dev.")
sumtab_ward_x <- xtable(sumtab_ward, caption = "Summary Statistics for Ward-Level Analysis (Gauteng Province)",
                        type = "latex")
digits(sumtab_ward_x) <- c(rep(0, 2), rep(2, 6))
print(sumtab_ward_x,
      caption.placement = 'top',
      scalebox = 1, 
      sanitize.text.function = function(x){x},
      label = "tab:wardsummarystats",
      file = "wardsummarystats.tex"
      )

## Appendix Table 2: Summary Statistics for Municipal-Level Analysis ##

# Councillor-level variables
dcoun <- d[, c("renom", "newprom", "years_incumbent", "ever_switchedparty")]
sumtab_municoun <- t(stat.desc(dcoun))[, c("nbr.val", "min", "max", "range", "median",
                                           "mean", "std.dev")]
varnames <- c("Renomination", "Promotion", "Years Incumbent", "Switched Party")
row.names(sumtab_municoun) <- varnames
# Municipal-level variables
dmuni_means <- ddply(d, .(cat_b), summarize,
                     service_pct_change_index1_MUN = mean(service_pct_change_index1),
                     service_rating_index2_MUN = mean(service_rating_index2),
                     AuditOpinion_MUN = mean(AuditOpinion),
                     formaldwell_pct_delta_MUN = mean(formaldwell_pct_delta),
                     logpopchange_MUN = mean(logpopchange),
                     black_avg_MUN = mean(black_avg),
                     post_secondary_avg_MUN = mean(post_secondary_avg),
                     win_margin_MUN = mean (win_margin))
sumtab_muni <- t(stat.desc(dmuni_means[, - 1]))[, c("nbr.val", "min", "max", "range", "median",
                                                    "mean", "std.dev")]
varnames <- c("Service Coverage $\\Delta$",
              "Service Rating Index",
              "Audit Opinion",
              "Formal Housing $\\Delta$",
              "Post Secondary (\\%)",
              "Ethnicity - African (\\%)",
              "Log Population $\\Delta$",
              "Win Margin (2011)")
rownames(sumtab_muni) <- varnames
sumtab_muni_comb <- rbind(sumtab_municoun, sumtab_muni)
colnames(sumtab_muni_comb) <- c("N obs", "Min", "Max", "Range", "Median", "Mean", "Std. Dev.")
sumtab_muni_comb_x <- xtable(sumtab_muni_comb, caption = "Summary Statistics for Municipal-Level Analysis (National)", type = "latex")
digits(sumtab_muni_comb_x) <- c(rep(0, 2), rep(2, 6))
print(sumtab_muni_comb_x,
      caption.placement = 'top',
      scalebox = 1,
      floating = TRUE,
      sanitize.text.function = function(x){x},
      label = "tab:munsummarystatsupdated",
      file = "munsummarystats_updated.tex")

## Appendix Table 3: Citizen Satisfaction, Councillor Renomination, and Political Competition in Gauteng (Logit) ##

# Column 1: All wards, continuous vote margin
all.modx <- glm(renom ~
                  local_satisfied_15*win_margin2011 
                + subj_sp_index15
                + obj_sp_delta_index_1315
                + formal_housing_delta_1315
                + unemp_delta_1315_reverse
                + high_income_avg
                + post_secondary_avg
                + black_avg
                + ce_index15
                + logpop2011
                + years_incumbent
                + ever_switchedparty 
                + I(Municipality),
                data = gcro_noBE, family = binomial(link = "logit"))
# Column 2: ANC wards, continuous vote margin
anc.modx <- glm(renom ~
                  local_satisfied_15*win_margin2011
                + subj_sp_index15
                + obj_sp_delta_index_1315
                + formal_housing_delta_1315
                + unemp_delta_1315_reverse
                + high_income_avg
                + post_secondary_avg
                + black_avg
                + ce_index15
                + logpop2011
                + years_incumbent
                + ever_switchedparty 
                + I(Municipality),
                data = gcro_ANC, family = binomial(link = "logit"))
# Column 3: DA wards, continuous vote margin
da.modx <- glm(renom ~
                 local_satisfied_15*win_margin2011
               + subj_sp_index15
               + obj_sp_delta_index_1315
               + formal_housing_delta_1315
               + unemp_delta_1315_reverse
               + subj_sp_index15 
               + high_income_avg
               + post_secondary_avg
               + black_avg
               + ce_index15
               + logpop2011
               + years_incumbent
               + ever_switchedparty 
               + I(Municipality),
               data = gcro_DA, family = binomial(link = "logit"))
# Create Table
labs <- c("Councillor Satis. x Win Margin", 
          "Councillor Satisfaction", "Win Margin (2011)",
          "Service Satisfaction Index",  "Service Coverage $\\Delta$", "Formal Housing $\\Delta$", "Employment $\\Delta$",
          "High Income", "Post Secondary Education", "Ethnicity - African", "Civic Engagement",
          "Log Population (2011)", "Years Incumbent", "Switched Party")
stargazer(all.modx,
          anc.modx,
          da.modx,
          covariate.labels = labs,
          column.sep.width = "1pt",
          no.space = TRUE,
          font.size = "footnotesize",
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes = "*$p<0.05$",
          notes.label = "",
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = F,
          column.labels = c("All", "ANC", "DA"),
          omit = c("I\\(", "Constant"),
          order = c(23, 1:13),
          omit.stat = c("ll", "aic"),
          add.lines = list(c("Municipality FE", "Yes", "Yes", "Yes")),
          label = "tab:wardcompetition",
          title = "Citizen Satisfaction, Councillor Renomination, and Political Competition in Gauteng (Logit)",
          type = "latex"
          ,
          out = "wardcompetition.tex"
          )

## Appendix Table 4: Councillor Renomination and Political Competition - All Municipalities (Logit) ##

# Column 1: Service rating and win margin interaction - all ruling
munmod.allx <- glm(renom ~
                     service_rating_index2*win_margin 
                   + service_pct_change_index1 
                   + formaldwell_pct_delta  
                   + logpopchange
                   + black_avg
                   + post_secondary_avg
                   + Ward
                   + years_incumbent
                   + ever_switchedparty 
                   + I(Province)
                   ,
                   data = drule, family = binomial(link = "logit"))
munmod.allx.c <- coeftest(munmod.allx, vcov = cluster.vcov(munmod.allx, cluster = drule$cat_b))
# Column 2: Service rating and win margin interaction - ANC only
munmod.ancx <- glm(renom ~
                     service_rating_index2*win_margin 
                   + service_pct_change_index1 
                   + formaldwell_pct_delta  
                   + logpopchange
                   + black_avg
                   + post_secondary_avg
                   + Ward
                   + years_incumbent
                   + ever_switchedparty 
                   + I(Province)
                   , 
                   data = drule_ANC, family = binomial(link = "logit"))
munmod.ancx.c <- coeftest(munmod.ancx, vcov = cluster.vcov(munmod.ancx, cluster = drule_ANC$cat_b))
# Column 3: Service rating and win margin interaction - non-ruling
munmod.nonrulex <- glm(renom ~
                         service_rating_index2*win_margin 
                       + service_pct_change_index1  
                       + formaldwell_pct_delta  
                       + logpopchange
                       + black_avg
                       + post_secondary_avg
                       + Ward
                       + years_incumbent
                       + ever_switchedparty 
                       + I(Province)
                       , 
                       data = dnonrule, family = binomial(link = "logit"))
munmod.nonrulex.c <- coeftest(munmod.nonrulex, vcov = cluster.vcov(munmod.nonrulex, cluster = dnonrule$cat_b))
# Create table
labs <- c("Svc. Rating Index x Win Margin", 
          "Service Rating Index", "Win Margin (2011)", "Service Coverage $\\Delta$",
          "Formal Housing $\\Delta$", "Log Population $\\Delta$",
          "Ethnicity - African (\\%)", "Post Secondary (\\%)", "Ward Councillor",
          "Years Incumbent", "Switched Party")          
ncounc <- c(nobs(munmod.allx), nobs(munmod.ancx), nobs(munmod.nonrulex))
nmunis <- length(unique(drule$cat_b))
nmunis_anc <- length(unique(drule_ANC$cat_b))
stargazer(munmod.allx.c, munmod.ancx.c, munmod.nonrulex.c,
          no.space = TRUE,
          column.sep.width = "0.1pt",
          covariate.labels = labs,
          column.labels = c("All Ruling", "ANC", "Non-Ruling"),
          omit = c("Province", "Constant"),
          order = c(19, c(1:10)), 
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = FALSE,
          font.size = "footnotesize",
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes = "*$p<0.05$",
          add.lines = list(
            c("N Councillors", ncounc),
            c("N Municipalities", rep(c(nmunis, nmunis_anc, nmunis), 2)),
            c("Province FE", rep("Yes", 6))
          ),
          notes.label = "RSE clustered at municipal level",
          title = "Councillor Renomination and Political Competition - All Municipalities (Logit)",
          label = "tab:muncompetitionupdated"
         ,
         out = "muncompetition_updated.tex"
        )

## Appendix Table 5: Determinants of Councillor Renomination - Top 10 PR Councillors, All Municipalities (Logit) ##

# Column 1: Top 10 party councillors, subjective ratings only
munmod.top10 <- glm(renom ~ 
                      service_rating_index2
                    + logpopchange
                    + black_avg
                    + post_secondary_avg
                    + win_margin
                    + years_incumbent
                    + ever_switchedparty 
                    + I(Province),
                    data = drule_top10, family = binomial(link = "logit"))
munmod.top10.c <- coeftest(munmod.top10, vcov = cluster.vcov(munmod.top10, drule_top10$cat_b))
# Column 2: Top 10 party councillors, ANC only, subjective ratings only
munmod.anc.top10 <- glm(renom ~ 
                          service_rating_index2
                        + logpopchange
                        + black_avg
                        + post_secondary_avg
                        + win_margin
                        + years_incumbent
                        + ever_switchedparty 
                        + I(Province),
                        data = drule_ANC_top10, family = binomial(link = "logit"))
munmod.anc.top10.c <- coeftest(munmod.anc.top10, vcov = cluster.vcov(munmod.anc.top10, drule_ANC_top10$cat_b))
# Column 3: Top 10 party councillors, non-ruling parties only, subjective ratings only (placebo)
munmod.nonrule.top10 <- glm(renom ~ 
                              service_rating_index2
                            + logpopchange
                            + black_avg
                            + post_secondary_avg
                            + win_margin
                            + years_incumbent
                            + ever_switchedparty 
                            + I(Province),
                            data = dnonrule_top10, family = binomial(link = "logit"))
munmod.nonrule.top10.c <- coeftest(munmod.nonrule.top10, vcov = cluster.vcov(munmod.nonrule.top10, dnonrule_top10$cat_b))
# Column 4: Top 10 party councillors, subjective and objective
munmod.obj.top10 <- glm(renom ~ 
                          service_rating_index2
                        + service_pct_change_index1
                        + formaldwell_pct_delta
                        + logpopchange
                        + black_avg
                        + post_secondary_avg
                        + win_margin
                        + years_incumbent
                        + ever_switchedparty 
                        + I(Province),
                        data = drule_top10, family = binomial(link = "logit"))
munmod.obj.top10.c <- coeftest(munmod.obj.top10, vcov = cluster.vcov(munmod.obj.top10, drule_top10$cat_b))
# Column 5: Top 10 party councillors, ANC only, subjective ratings only
munmod.anc.obj.top10 <- glm(renom ~ 
                              service_rating_index2
                            + service_pct_change_index1
                            + formaldwell_pct_delta
                            + logpopchange
                            + black_avg
                            + post_secondary_avg
                            + win_margin
                            + years_incumbent
                            + ever_switchedparty 
                            + I(Province),
                            data = drule_ANC_top10, family = binomial(link = "logit"))
munmod.anc.obj.top10.c <- coeftest(munmod.anc.obj.top10, vcov = cluster.vcov(munmod.anc.obj.top10, drule_ANC_top10$cat_b))
# Column 6: Top 10 party councillors, non-ruling parties only, subjective ratings only (placebo)
munmod.nonrule.obj.top10 <- glm(renom ~ 
                                  service_rating_index2
                                + service_pct_change_index1
                                + formaldwell_pct_delta
                                + logpopchange
                                + black_avg
                                + post_secondary_avg
                                + win_margin
                                + years_incumbent
                                + ever_switchedparty 
                                + I(Province),
                                data = dnonrule_top10, family = binomial(link = "logit"))
munmod.nonrule.obj.top10.c <- coeftest(munmod.nonrule.obj.top10, vcov = cluster.vcov(munmod.nonrule.obj.top10, dnonrule_top10$cat_b))

# Create table
labs <- c("Service Rating Index", "Service Coverage $\\Delta$", 
          "Formal Housing $\\Delta$", "Log Population $\\Delta$",
          "Ethnicity - African (\\%)", "Post Secondary (\\%)", "Win Margin (2011)",
          "Years Incumbent", "Switched Party")
nmunis <- length(unique(drule_top10$cat_b))
nmunis_anc <- length(unique(drule_ANC_top10$cat_b))
stargazer(munmod.top10.c,
          munmod.anc.top10.c,
          munmod.nonrule.top10.c,
          munmod.obj.top10.c,
          munmod.anc.obj.top10.c,
          munmod.nonrule.obj.top10.c,
          no.space = TRUE,
          covariate.labels = labs,
          column.labels = rep(c("All Ruling", "ANC", "Non-Ruling"), 2),
          omit = c("Province", "Constant"),
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = FALSE,
          add.lines = list(
            c("N Councillors", nobs(munmod.top10), nobs(munmod.anc.top10), nobs(munmod.nonrule.top10), nobs(munmod.obj.top10),
              nobs(munmod.obj.top10), nobs(munmod.nonrule.obj.top10)),
            c("N Municipalities", rep(c(nmunis, nmunis_anc, nmunis), 2)),
            c("Province FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")
          ),
          font.size = "footnotesize",
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes = "*$p < 0.05$",
          notes.label = "RSE clustered at municipal level",
          title = "Determinants of Councillor Renomination - Top 10 PR Councillors, All Municipalities (Logit)",
          label = "tab:munmaintop10"
          ,
          out = "munmain.tex"
          )

# Substantive coefficient
munres_top10rule <- mainmunFun(munmod.obj.top10, drule_top10)
munres_top10rule

## Appendix Table 6: Councillor Renomination and Political Competition - Top 10 PR Councillors, All Municipalities (Logit) ##

# Column 1: Top 10 party councillors, subjective and objective
munmodx.obj.top10 <- glm(renom ~ 
                           service_rating_index2*win_margin
                         + service_pct_change_index1
                         + formaldwell_pct_delta
                         + logpopchange
                         + black_avg
                         + post_secondary_avg
                         + win_margin
                         + Ward
                         + years_incumbent
                         + ever_switchedparty 
                         + I(Province),
                         data = drule_top10, family = binomial(link = "logit"))
munmodx.obj.top10.c <- coeftest(munmodx.obj.top10, vcov = cluster.vcov(munmodx.obj.top10, drule_top10$cat_b))
# Column 2: Top 10 party councillors, ANC only
munmodx.anc.obj.top10 <- glm(renom ~ 
                               service_rating_index2*win_margin
                             + service_pct_change_index1
                             + formaldwell_pct_delta
                             + logpopchange
                             + black_avg
                             + post_secondary_avg
                             + win_margin
                             + Ward
                             + years_incumbent
                             + ever_switchedparty 
                             + I(Province),
                             data = drule_ANC_top10, family = binomial(link = "logit"))
munmodx.anc.obj.top10.c <- coeftest(munmodx.anc.obj.top10, vcov = cluster.vcov(munmodx.anc.obj.top10, drule_ANC_top10$cat_b))
# Column 3: Top 10 party councillors, non-ruling parties
munmodx.nonrule.obj.top10 <- glm(renom ~ 
                                   service_rating_index2*win_margin
                                 + service_pct_change_index1
                                 + formaldwell_pct_delta
                                 + logpopchange
                                 + black_avg
                                 + post_secondary_avg
                                 + win_margin
                                 + Ward
                                 + years_incumbent
                                 + ever_switchedparty 
                                 + I(Province),
                                 data = dnonrule_top10, family = binomial(link = "logit"))
munmodx.nonrule.obj.top10.c <- coeftest(munmodx.nonrule.obj.top10, vcov = cluster.vcov(munmodx.nonrule.obj.top10, dnonrule_top10$cat_b))

# Create table
labs <- c("Svc. Rating Index x Win Margin", 
          "Service Rating Index", "Win Margin (2011)", "Service Coverage $\\Delta$",
          "Formal Housing $\\Delta$", "Log Population $\\Delta$",
          "Ethnicity - African (\\%)", "Post Secondary (\\%)", "Ward Councillor",
          "Years Incumbent", "Switched Party")          
ncounc <- c(nobs(munmodx.obj.top10), nobs(munmodx.anc.obj.top10), nobs(munmodx.nonrule.obj.top10))
nmunis <- length(unique(drule_top10$cat_b))
nmunis_anc <- length(unique(drule_ANC_top10$cat_b))
stargazer(munmodx.obj.top10.c, munmodx.anc.obj.top10.c, munmodx.nonrule.obj.top10.c,
          no.space = TRUE,
          column.sep.width = "0.1pt",
          covariate.labels = labs,
          column.labels = c("All Ruling", "ANC", "Non-Ruling"),
          omit = c("Province", "Constant"),
          order = c(18, c(1:9)), 
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = FALSE,
          font.size = "footnotesize",
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes = "*$p<0.05$",
          add.lines = list(
            c("N Councillors", ncounc),
            c("N Municipalities", rep(c(nmunis, nmunis_anc, nmunis), 2)),
            c("Province FE", rep("Yes", 6))
          ),
          notes.label = "RSE clustered at municipal level",
          title = "Councillor Renomination and Political Competition - Top 10 PR Councillors, All Municipalities (Logit)",
          label = "tab:muncompetitiontop10"
          ,
          out = "muncompetition_top10.tex"
)

## Appendix Table 7: Determinants of Councillor Renomination in Gauteng (Logit) - ANC-ruled Municipalities Only ##

# Column 1: All Wards - satisfaction only
all.mod.noMV <- glm(renom ~
                      local_satisfied_15 
                    + high_income_avg
                    + post_secondary_avg
                    + black_avg
                    + ce_index15
                    + logpop2011
                    + win_margin2011
                    + years_incumbent
                    + ever_switchedparty 
                    + I(Municipality),
                    data = gcro_noMV, family = binomial(link = "logit"))
# Column 2: All Wards - councillor satisfaction and service satisfaction
all.mod.subj.noMV <- glm(renom ~ 
                           local_satisfied_15
                         + subj_sp_index15
                         + high_income_avg
                         + post_secondary_avg
                         + black_avg
                         + ce_index15
                         + logpop2011
                         + win_margin2011
                         + years_incumbent
                         + ever_switchedparty 
                         + I(Municipality),
                         data = gcro_noMV, family = binomial(link = "logit"))
# Column 3: All wards - councillor satisfaction, service satisfaction, objective services
all.mod.obj.noMV <- glm(renom ~ 
                          local_satisfied_15 
                        + subj_sp_index15
                        + obj_sp_delta_index_1315
                        + formal_housing_delta_1315
                        + unemp_delta_1315_reverse
                        + high_income_avg
                        + post_secondary_avg
                        + black_avg
                        + ce_index15
                        + logpop2011
                        + win_margin2011
                        + years_incumbent
                        + ever_switchedparty 
                        + I(Municipality),
                        data = gcro_noMV, family = binomial(link = "logit"))
# Column 4: Councillor satisfaction only, ANC wards only
anc.mod.noMV <- glm(renom ~ 
                      local_satisfied_15
                    + high_income_avg
                    + post_secondary_avg
                    + black_avg
                    + ce_index15
                    + logpop2011
                    + win_margin2011
                    + years_incumbent
                    + ever_switchedparty 
                    + I(Municipality),
                    data = gcro_ANC_muns, family = binomial(link = "logit"))
# Column 5: councillor satisfaction and service satisfaction, ANC only
anc.mod.subj.noMV <- glm(renom ~
                           local_satisfied_15
                         + subj_sp_index15
                         + high_income_avg
                         + post_secondary_avg
                         + black_avg
                         + ce_index15
                         + logpop2011
                         + win_margin2011
                         + years_incumbent
                         + ever_switchedparty 
                         + I(Municipality),
                         data = gcro_ANC_muns, family = binomial(link = "logit"))
# Column 6: councillor satisfaction, service satisfaction, and objective sp, ANC only
anc.mod.obj.noMV <- glm(renom ~ 
                          local_satisfied_15
                        + subj_sp_index15
                        + obj_sp_delta_index_1315
                        + formal_housing_delta_1315
                        + unemp_delta_1315_reverse
                        + high_income_avg
                        + post_secondary_avg
                        + black_avg
                        + ce_index15
                        + logpop2011
                        + win_margin2011
                        + years_incumbent
                        + ever_switchedparty 
                        + I(Municipality),
                        data = gcro_ANC_muns, family = binomial(link = "logit"))
# Column 7: councillor satisfaction only, DA wards
da.mod.noMV <- glm(renom ~ 
                     local_satisfied_15
                   + high_income_avg
                   + post_secondary_avg
                   + black_avg
                   + ce_index15
                   + logpop2011
                   + win_margin2011
                   + years_incumbent
                   + ever_switchedparty 
                   + I(Municipality),
                   data = gcro_DA_noMV, family = binomial(link = "logit"))
# Column 8: councillor satisfaction and service satisfaction, DA only
da.mod.subj.noMV <- glm(renom ~ 
                          local_satisfied_15
                        + subj_sp_index15
                        + high_income_avg
                        + post_secondary_avg
                        + black_avg
                        + ce_index15
                        + logpop2011
                        + win_margin2011
                        + years_incumbent
                        + ever_switchedparty 
                        + I(Municipality),
                        data = gcro_DA_noMV, family = binomial(link = "logit"))
# Column 9: councillor satisfaction, service satisfaction, and objective SP, DA only
da.mod.obj.noMV <- glm(renom ~
                         local_satisfied_15
                       + subj_sp_index15
                       + obj_sp_delta_index_1315
                       + formal_housing_delta_1315
                       + unemp_delta_1315_reverse
                       + high_income_avg
                       + post_secondary_avg
                       + black_avg
                       + ce_index15
                       + logpop2011
                       + win_margin2011
                       + years_incumbent
                       + ever_switchedparty 
                       + I(Municipality),
                       data = gcro_DA_noMV, family = binomial(link = "logit"))
# Create table
labs <- c("Councillor Satisfaction", "Service Satisfaction Index",
          "Service Coverage $\\Delta$", "Formal Housing $\\Delta$",
          "Employment $\\Delta$", "High Income (\\%)", "Post Secondary (\\%)",
          "Ethnicity - African (\\%)", "Civic Engagement",
          "Log Population (2011)", "Win Margin (2011)", "Years Incumbent", "Switched Party")
stargazer(all.mod.noMV,
          all.mod.subj.noMV,
          all.mod.obj.noMV,
          anc.mod.noMV,
          anc.mod.subj.noMV,
          anc.mod.obj.noMV,
          da.mod.noMV,
          da.mod.subj.noMV,
          da.mod.obj.noMV,
          omit = c("I\\(", "Constant"),
          no.space = TRUE,
          covariate.labels = labs,
          column.labels = c(rep("All", 3), rep("ANC", 3), rep("DA", 3)),
          omit.stat = c("ll", "aic"),
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = F,
          font.size = "footnotesize",
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes.label = "",
          notes = "*$p<0.05$",
          title = "Determinants of Councillor Renomination in Gauteng (Logit) - ANC-ruled Municipalities Only",
          add.lines = list(c("Municipality FE", rep("Yes", 9))),
          out = "wardmainnoMV.tex",
          label = "tab:wardmainnoMV"
)

## Appendix Table 8: Determinants of Councillor Renomination in Gauteng (Logit) - Including Protest Participation ##

# Column 1: All Wards - satisfaction only
all.mod.protest <- glm(renom ~ 
                         local_satisfied_15 
                       + high_income_avg
                       + post_secondary_avg
                       + black_avg
                       + participate_protest_15
                       + logpop2011
                       + win_margin2011
                       + years_incumbent
                       + ever_switchedparty 
                       + I(Municipality),
                       data = gcro_noBE, family = binomial(link = "logit"))
# Column 2: All Wards - councillor satisfaction and service satisfaction
all.mod.subj.protest <- glm(renom ~ 
                              local_satisfied_15 
                            + subj_sp_index15
                            + high_income_avg
                            + post_secondary_avg
                            + black_avg
                            + participate_protest_15
                            + logpop2011
                            + win_margin2011
                            + years_incumbent
                            + ever_switchedparty 
                            + I(Municipality),
                            data = gcro_noBE, family = binomial(link = "logit"))
# Column 3: All wards - councillor satisfaction, service satisfaction, objective services
all.mod.obj.protest <- glm(renom ~
                             local_satisfied_15 
                           + subj_sp_index15
                           + obj_sp_delta_index_1315
                           + formal_housing_delta_1315
                           + unemp_delta_1315_reverse
                           + high_income_avg
                           + post_secondary_avg
                           + black_avg
                           + participate_protest_15
                           + logpop2011
                           + win_margin2011
                           + years_incumbent
                           + ever_switchedparty 
                           + I(Municipality),
                           data = gcro_noBE, family = binomial(link = "logit"))
# Column 4: Councillor satisfaction only, ANC wards only
anc.mod.protest <- glm(renom ~ 
                         local_satisfied_15
                       + high_income_avg
                       + post_secondary_avg
                       + black_avg
                       + participate_protest_15
                       + logpop2011
                       + win_margin2011
                       + years_incumbent
                       + ever_switchedparty 
                       + I(Municipality),
                       data = gcro_ANC, family = binomial(link = "logit"))
# Column 5: councillor satisfaction and service satisfaction, ANC only
anc.mod.subj.protest <- glm(renom ~ 
                              local_satisfied_15
                            + subj_sp_index15
                            + high_income_avg
                            + post_secondary_avg
                            + black_avg
                            + participate_protest_15
                            + logpop2011
                            + win_margin2011
                            + years_incumbent
                            + ever_switchedparty 
                            + I(Municipality),
                            data = gcro_ANC, family = binomial(link = "logit"))
# Column 6: councillor satisfaction, service satisfaction, and objective sp, ANC only
anc.mod.obj.protest <- glm(renom ~ 
                             local_satisfied_15
                           + subj_sp_index15
                           + obj_sp_delta_index_1315
                           + formal_housing_delta_1315
                           + unemp_delta_1315_reverse
                           + high_income_avg
                           + post_secondary_avg
                           + black_avg
                           + participate_protest_15
                           + logpop2011
                           + win_margin2011
                           + years_incumbent
                           + ever_switchedparty 
                           + I(Municipality),
                           data = gcro_ANC, family = binomial(link = "logit"))
# Column 7: councillor satisfaction only, DA wards
da.mod.protest <- glm(renom ~
                        local_satisfied_15
                      + high_income_avg
                      + post_secondary_avg
                      + black_avg
                      + participate_protest_15
                      + logpop2011
                      + win_margin2011
                      + years_incumbent
                      + ever_switchedparty 
                      + I(Municipality),
                      data = gcro_DA, family = binomial(link = "logit"))
# Column 8: councillor satisfaction and service satisfaction, DA only
da.mod.subj.protest <- glm(renom ~
                             local_satisfied_15
                           + subj_sp_index15
                           + high_income_avg
                           + post_secondary_avg
                           + black_avg
                           + participate_protest_15
                           + logpop2011
                           + win_margin2011
                           + years_incumbent
                           + ever_switchedparty 
                           + I(Municipality),
                           data = gcro_DA, family = binomial(link = "logit"))
summary(da.mod.subj.protest)
# Column 9: councillor satisfaction, service satisfaction, and objective SP, DA only
da.mod.obj.protest <- glm(renom ~ 
                            local_satisfied_15
                          + subj_sp_index15
                          + obj_sp_delta_index_1315
                          + formal_housing_delta_1315
                          + unemp_delta_1315_reverse
                          + high_income_avg
                          + post_secondary_avg
                          + black_avg
                          + participate_protest_15
                          + logpop2011
                          + win_margin2011
                          + years_incumbent
                          + ever_switchedparty 
                          + I(Municipality),
                          data = gcro_DA, family = binomial(link = "logit"))
summary(da.mod.obj.protest)
# Create table
labs <- c("Councillor Satisfaction", "Service Satisfaction Index",
          "Service Coverage $\\Delta$", "Formal Housing $\\Delta$",
          "Employment $\\Delta$", "High Income (\\%)", "Post Secondary (\\%)",
          "Ethnicity - African (\\%)", "Protest Participation",
          "Log Population (2011)", "Win Margin (2011)", "Years Incumbent", "Switched Party")
stargazer(all.mod.protest,
          all.mod.subj.protest,
          all.mod.obj.protest,
          anc.mod.protest,
          anc.mod.subj.protest,
          anc.mod.obj.protest,
          da.mod.protest,
          da.mod.subj.protest,
          da.mod.obj.protest,
          omit = c("I\\(", "Constant"),
          no.space = TRUE,
          covariate.labels = labs,
          column.labels = c(rep("All", 3), rep("ANC", 3), rep("DA", 3)),
          omit.stat = c("ll", "aic"),
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = F,
          font.size = "footnotesize",
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes.label = "",
          notes = "*$p<0.05$",
          title = "Determinants of Councillor Renomination in Gauteng (Logit) - Including Protest Participation",
          add.lines = list(c("Municipality FE", rep("Yes", 9))),
          out = "wardmainprotest.tex",
          label = "tab:wardmainprotest"
)

## Appendix Table 9: Service Improvements and Councillor Renomination in Gauteng (Logit) ##

# Column 1: All wards - objective measures only
all.mod.objonly <- glm(renom ~ 
                        obj_sp_delta_index_1315
                       + formal_housing_delta_1315
                       + unemp_delta_1315_reverse
                       + high_income_avg
                       + post_secondary_avg
                       + black_avg
                       + ce_index15
                       + logpop2011
                       + win_margin2011
                       + years_incumbent
                       + ever_switchedparty 
                       + I(Municipality),
                       data = gcro_noBE, family = binomial(link = "logit"))
# Column 2: ANC wards only - objective measures only
anc.mod.objonly <- glm(renom ~ 
                         obj_sp_delta_index_1315
                       + formal_housing_delta_1315
                       + unemp_delta_1315_reverse
                       + high_income_avg
                       + post_secondary_avg
                       + black_avg
                       + ce_index15
                       + logpop2011
                       + win_margin2011
                       + years_incumbent
                       + ever_switchedparty 
                       + I(Municipality),
                       data = gcro_ANC_muns, family = binomial(link = "logit"))
# Column 3: DA wards only - objective measures only
da.mod.objonly <- glm(renom ~ 
                        obj_sp_delta_index_1315
                      + formal_housing_delta_1315
                      + unemp_delta_1315_reverse
                      + high_income_avg
                      + post_secondary_avg
                      + black_avg
                      + ce_index15
                      + logpop2011
                      + win_margin2011
                      + years_incumbent
                      + ever_switchedparty 
                      + I(Municipality),
                      data = gcro_DA_noMV, family = binomial(link = "logit"))
# Column 4: All wards, objective service interaction
osponly.modx <- glm(renom ~
                      obj_sp_delta_index_1315*win_margin2011
                    + formal_housing_delta_1315
                    + unemp_delta_1315_reverse
                    + high_income_avg
                    + post_secondary_avg
                    + black_avg
                    + ce_index15
                    + logpop2011
                    + years_incumbent
                    + ever_switchedparty 
                    + I(Municipality),
                    data = gcro_noBE, family = binomial(link = "logit"))
# Column 5: ANC wards only, objective service interaction
osponly.modx.anc <- glm(renom ~
                          obj_sp_delta_index_1315*win_margin2011
                        + formal_housing_delta_1315
                        + unemp_delta_1315_reverse
                        + high_income_avg
                        + post_secondary_avg
                        + black_avg
                        + ce_index15
                        + logpop2011
                        + years_incumbent
                        + ever_switchedparty 
                        + I(Municipality),
                        data = gcro_ANC_muns, family = binomial(link = "logit"))
# Column 6: DA wards only, objective service interaction
osponly.modx.da <- glm(renom ~
                         obj_sp_delta_index_1315*win_margin2011
                       + formal_housing_delta_1315
                       + unemp_delta_1315_reverse
                       + high_income_avg
                       + post_secondary_avg
                       + black_avg
                       + ce_index15
                       + logpop2011
                       + years_incumbent
                       + ever_switchedparty 
                       + I(Municipality),
                       data = gcro_DA_noMV, family = binomial(link = "logit"))
# Create table
labs <- c("Service Coverage $\\Delta$ x Win Margin",
          "Service Coverage $\\Delta$", "Win Margin", "Formal Housing $\\Delta$",
          "Employment $\\Delta$", "High Income (\\%)", "Post Secondary (\\%)",
          "Ethnicity - African (\\%)", "Civic Engagement",
          "Log Population (2011)", "Years Incumbent", "Switched Party")
stargazer(all.mod.objonly, anc.mod.objonly, da.mod.objonly,
          osponly.modx, osponly.modx.anc, osponly.modx.da,
          no.space = TRUE,
          covariate.labels = labs,
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes.label = "",
          notes = "*$p<0.05$",
          order = c(21, 1, 9, c(2:8), c(10:20)),
          column.labels = rep(c("All", "ANC", "DA"), 2),
          omit = c(" - ", "Constant"),
          omit.stat = c("ll", "aic"),
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = FALSE,
          title = "Service Improvement and Councillor Renomination in Gauteng (Logit)",
          add.lines = list(c("Municipality FE", rep("Yes", 6))),
          label = "tab:wardobj"
          ,
          out = "wardobj.tex")

## Appendix Table 10: Service Improvements and Councillor Renomination - All Municipalities (Logit) ##

# Column 1: All ruling party councillors, objective changes only
munmod.all.objonly <- glm(renom ~
                            service_pct_change_index1
                          + formaldwell_pct_delta  
                          + logpopchange
                          + black_avg
                          + post_secondary_avg
                          + win_margin
                          + Ward
                          + years_incumbent
                          + ever_switchedparty
                          + I(Province), 
                          data = drule, family = binomial(link = "logit"))
munmod.all.objonly.c <- coeftest(munmod.all.objonly, vcov = cluster.vcov(munmod.all.objonly, drule$cat_b))
# Column 2: ANC only, objective changes only
munmod.alla.objonly <- glm(renom ~
                             service_pct_change_index1
                           + formaldwell_pct_delta  
                           + logpopchange
                           + black_avg
                           + post_secondary_avg
                           + win_margin
                           + Ward
                           + years_incumbent
                           + ever_switchedparty
                           + I(Province), 
                           data = drule_ANC, family = binomial(link = "logit"))
munmod.alla.objonly.c <- coeftest(munmod.alla.objonly, vcov = cluster.vcov(munmod.alla.objonly, drule_ANC$cat_b))
# Column 3: Non-ruling party councillors only, objective changes only
munmod.nonrule.objonly <- glm(renom ~ 
                                service_pct_change_index1
                              + formaldwell_pct_delta  
                              + logpopchange
                              + black_avg
                              + post_secondary_avg
                              + win_margin
                              + Ward
                              + years_incumbent
                              + ever_switchedparty
                              + I(Province), 
                              data = dnonrule, family = binomial(link = "logit"))
munmod.nonrule.objonly.c <- coeftest(munmod.nonrule.objonly, vcov = cluster.vcov(munmod.nonrule.objonly, dnonrule$cat_b))
# Column 4: All ruling party councillors, service change interacted with political competition, not controlling for service ratings
munmod.osp.objonly.x <- glm(renom ~
                              service_pct_change_index1*win_margin 
                            + formaldwell_pct_delta  
                            + logpopchange
                            + black_avg
                            + post_secondary_avg
                            + Ward
                            + years_incumbent
                            + ever_switchedparty
                            + I(Province), 
                            data = drule, family = binomial(link = "logit"))
munmod.osp.objonly.x.c <- coeftest(munmod.osp.objonly.x, vcov = cluster.vcov(munmod.osp.objonly.x, cluster = drule$cat_b))
# Column 5: ANC only, service change interacted with political competition, not controlling for service ratings
munmod.osp.objonly.ancx <- glm(renom ~
                                 service_pct_change_index1*win_margin 
                               + formaldwell_pct_delta  
                               + logpopchange
                               + black_avg
                               + post_secondary_avg
                               + Ward
                               + years_incumbent
                               + ever_switchedparty 
                               + I(Province), 
                               data = drule_ANC, family = binomial(link = "logit"))
munmod.osp.objonly.ancx.c <- coeftest(munmod.osp.objonly.ancx, vcov = cluster.vcov(munmod.osp.objonly.ancx, cluster = drule_ANC$cat_b))
# Column 6: Non-ruling parties only, Service change interacted with political competition, not controlling for service ratings
munmod.osp.objonly.nonrulex <- glm(renom ~
                                     service_pct_change_index1*win_margin 
                                   + formaldwell_pct_delta  
                                   + logpopchange
                                   + black_avg
                                   + post_secondary_avg
                                   + Ward
                                   + years_incumbent
                                   + ever_switchedparty
                                   + I(Province), 
                                   data = dnonrule, family = binomial(link = "logit"))
munmod.osp.objonly.nonrulex.c <- coeftest(munmod.osp.objonly.nonrulex, vcov = cluster.vcov(munmod.osp.objonly.nonrulex, cluster = dnonrule$cat_b))
# Create table
labs <- c("Svc. Coverage $\\Delta$ x Win Margin", "Svc. Coverage $\\Delta$", "Win Margin",
          "Formal Housing $\\Delta$", "Log Population $\\Delta$",
          "Ethnicity - African (\\%)", "Post Secondary (\\%)", "Ward Councillor",
          "Years Incumbent", "Switched Party")
nmunis <- length(unique(drule$cat_b))
nmunis_anc <- length(unique(drule_ANC$cat_b))
stargazer(munmod.all.objonly.c, munmod.alla.objonly.c, munmod.nonrule.objonly.c,
          munmod.osp.objonly.x.c, munmod.osp.objonly.ancx.c, munmod.osp.objonly.nonrulex.c,
          no.space = TRUE,
          covariate.labels = labs,
          column.labels = rep(c("All Ruling", "ANC", "Non-Ruling"), 2),
          order = c(18, 1, 6, c(2:5)),
          omit = c("Province", "Constant"),
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes = "*$p<0.05$",
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = FALSE,
          add.lines = list(
            c("N Councillors", nobs(munmod.all.objonly), nobs(munmod.alla.objonly), nobs(munmod.nonrule.objonly),
              nobs(munmod.osp.objonly.x), nobs(munmod.osp.objonly.ancx), nobs(munmod.osp.objonly.nonrulex)),
            c("N Municipalities", rep(c(nmunis, nmunis_anc, nmunis), 2)),
            c("Province FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")
          ),
          notes.label = "RSE clustered at municipal level",
          title = "Service Improvement and Councillor Renomination - All Municipalities (Logit)",
          label = "tab:munobjupdated"
          ,
          out = "munobj_updated.tex"
          )

## Appendix Table 11: Determinants of Councillor Renomination - All Municipalities (Logit) without Province FE ## 

# Column 1: All ruling party councillors, subjective ratings only
munmod.all.nofe <- glm(renom ~
                         service_rating_index2
                       + logpopchange
                       + black_avg
                       + post_secondary_avg
                       + win_margin
                       + Ward
                       + years_incumbent
                       + ever_switchedparty,
                  data = drule, family = binomial(link = "logit"))
munmod.all.nofe.c <- coeftest(munmod.all.nofe, vcov = cluster.vcov(munmod.all.nofe, drule$cat_b))
# Column 2: ANC only, subjective ratings only
munmod.anc.nofe <- glm(renom ~
                         service_rating_index2
                       + logpopchange
                       + black_avg
                       + post_secondary_avg
                       + win_margin
                       + Ward
                       + years_incumbent
                       + ever_switchedparty,
                       data = drule_ANC, family = binomial(link = "logit"))
munmod.anc.nofe.c <- coeftest(munmod.anc.nofe, vcov = cluster.vcov(munmod.anc.nofe, drule_ANC$cat_b))
# Column 3: Non-ruling party councillors only, subjective ratings only
munmod.nonrule.nofe <- glm(renom ~
                             service_rating_index2
                           + logpopchange
                           + black_avg
                           + post_secondary_avg
                           + win_margin
                           + Ward
                           + years_incumbent
                           + ever_switchedparty, 
                      data = dnonrule, family = binomial(link = "logit"))
munmod.nonrule.nofe.c <- coeftest(munmod.nonrule.nofe, vcov = cluster.vcov(munmod.nonrule.nofe, dnonrule$cat_b))
# Column 4: All ruling party councillors, controlling for objective service improvements
munmod.all.obj.nofe <- glm(renom ~
                             service_rating_index2
                           + service_pct_change_index1
                           + formaldwell_pct_delta
                           + logpopchange
                           + black_avg
                           + post_secondary_avg
                           + win_margin
                           + Ward
                           + years_incumbent
                           + ever_switchedparty,
                           data = drule, family = binomial(link = "logit"))
munmod.all.obj.nofe.c <- coeftest(munmod.all.obj.nofe, vcov = cluster.vcov(munmod.all.obj.nofe, drule$cat_b))
# Column 5: ANC councillors only, controlling for objective service improvements
munmod.anc.obj.nofe <- glm(renom ~
                             service_rating_index2
                           + service_pct_change_index1
                           + formaldwell_pct_delta
                           + logpopchange
                           + black_avg
                           + post_secondary_avg
                           + win_margin
                           + Ward
                           + years_incumbent
                           + ever_switchedparty 
                       , data = drule_ANC, family = binomial(link = "logit"))
munmod.anc.obj.nofe.c <- coeftest(munmod.anc.obj.nofe, vcov = cluster.vcov(munmod.anc.obj.nofe, drule_ANC$cat_b))
# Column 6: Non-ruling party councillors only, controlling for objective service improvements
munmod.nonrule.obj.nofe <- glm(renom ~
                                 service_rating_index2
                               + service_pct_change_index1
                               + formaldwell_pct_delta
                               + logpopchange
                               + black_avg
                               + post_secondary_avg
                               + win_margin
                               + Ward
                               + years_incumbent
                               + ever_switchedparty,
                               data = dnonrule, family = binomial(link = "logit"))
munmod.nonrule.obj.nofe.c <- coeftest(munmod.nonrule.obj.nofe, vcov = cluster.vcov(munmod.nonrule.obj.nofe, dnonrule$cat_b))
# Create table
labs <- c("Service Rating Index", "Service Coverage $\\Delta$", 
          "Formal Housing $\\Delta$", "Log Population $\\Delta$",
          "Ethnicity - African (\\%)", "Post Secondary (\\%)", "Win Margin (2011)", "
          Ward Councillor",
          "Years Incumbent", "Switched Party")
nmunis <- length(unique(drule$cat_b))
nmunis_anc <- length(unique(drule_ANC$cat_b))
stargazer(munmod.all.nofe.c,
          munmod.anc.nofe.c,
          munmod.nonrule.nofe.c,
          munmod.all.obj.nofe.c,
          munmod.anc.obj.nofe.c,
          munmod.nonrule.obj.nofe.c,
          no.space = TRUE,
          covariate.labels = labs,
          column.labels = rep(c("All Ruling", "ANC", "Non-Ruling"), 2),
          omit = c("Province", "Constant"),
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = FALSE,
          add.lines = list(
            c("N Councillors", nobs(munmod.all.nofe), nobs(munmod.anc.nofe), nobs(munmod.nonrule.nofe), nobs(munmod.all.obj.nofe),
              nobs(munmod.anc.obj.nofe), nobs(munmod.nonrule.obj.nofe)),
            c("N Municipalities", rep(c(nmunis, nmunis_anc, nmunis), 2)),
            c("Province FE", "No", "No", "No", "No", "No", "No")
          ),
          font.size = "footnotesize",
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes = "*$p < 0.05$",
          notes.label = "RSE clustered at municipal level",
          title = "Determinants of Councillor Renomination - All Municipalities (Logit) without Province FE",
          label = "tab:munmainnofeupdated"
          ,
          out = "munmainnofe_updated.tex"
          )

## Appendix Table 12: Determinants of Councillor Renomination - All Municipalities (Logit), Ward vs. PR ##

# Column 1: All ruling party WARD councillors, controlling for objective service improvements
munmod.all.obj.ward <- glm(renom ~
                             service_rating_index2
                           + service_pct_change_index1
                           + formaldwell_pct_delta  
                           + logpopchange
                           + black_avg
                           + post_secondary_avg
                           + win_margin
                           + years_incumbent
                           + ever_switchedparty 
                           + I(Province)
                           , data = drule_ward, family = binomial(link = "logit"))
munmod.all.obj.ward.c <- coeftest(munmod.all.obj.ward, vcov = cluster.vcov(munmod.all.obj.ward, drule_ward$cat_b))
# Column 2: ANC WARD councillors only, controlling for objective service improvements
munmod.alla.obj.ward <- glm(renom ~  
                              service_rating_index2
                            + service_pct_change_index1 
                            + formaldwell_pct_delta  
                            + logpopchange 
                            + black_avg 
                            + post_secondary_avg 
                            + win_margin
                            + years_incumbent 
                            + ever_switchedparty 
                            + I(Province)
                            , data = drule_ANC[which(drule_ANC$Ward == 1), ], family = binomial(link = "logit"))
munmod.alla.obj.ward.c <- coeftest(munmod.alla.obj.ward, vcov = cluster.vcov(munmod.alla.obj.ward, drule_ANC$cat_b[which(drule_ANC$Ward == 1)]))
# Column 3: Non-ruling party WARD councillors only, controlling for objective service improvements
munmod.nonrule.obj.ward <- glm(renom ~ 
                                 service_rating_index2
                               + service_pct_change_index1 
                               + formaldwell_pct_delta  
                               + logpopchange 
                               + black_avg 
                               + post_secondary_avg 
                               + win_margin
                               + years_incumbent 
                               + ever_switchedparty 
                               + I(Province)
                               , data = dnonrule[which(dnonrule$Ward == 1), ], family = binomial(link = "logit"))
munmod.nonrule.obj.ward.c <- coeftest(munmod.nonrule.obj.ward, vcov = cluster.vcov(munmod.nonrule.obj.ward, dnonrule$cat_b[which(dnonrule$Ward ==1)]))
# Column 4: All ruling party PR councillors, controlling for objective service improvements
munmod.all.obj.pr <- glm(renom ~ 
                           service_rating_index2
                         + service_pct_change_index1
                         + formaldwell_pct_delta  
                         + logpopchange
                         + black_avg
                         + post_secondary_avg
                         + win_margin
                         + years_incumbent 
                         + ever_switchedparty 
                         + I(Province)
                         , data = drule_pr, family = binomial(link = "logit"))
munmod.all.obj.pr.c <- coeftest(munmod.all.obj.pr, vcov = cluster.vcov(munmod.all.obj.pr, drule_pr$cat_b))
# Column 5: ANC PR councillors only, controlling for objective service improvements
munmod.alla.obj.pr <- glm(renom ~  
                            service_rating_index2
                          + service_pct_change_index1
                          + formaldwell_pct_delta  
                          + logpopchange
                          + black_avg
                          + post_secondary_avg
                          + win_margin
                          + years_incumbent
                          + ever_switchedparty 
                          + I(Province)
                          , data = drule_ANC[which(drule_ANC$PR == 1), ], family = binomial(link = "logit"))
munmod.alla.obj.pr.c <- coeftest(munmod.alla.obj.pr, vcov = cluster.vcov(munmod.alla.obj.pr, drule_ANC$cat_b[which(drule_ANC$PR == 1)]))
# Column 6: Non-ruling party PR councillors only, controlling for objective service improvements
munmod.nonrule.obj.pr <- glm(renom ~ 
                               service_rating_index2
                             + service_pct_change_index1 
                             + formaldwell_pct_delta  
                             + logpopchange 
                             + black_avg 
                             + post_secondary_avg 
                             + win_margin
                             + years_incumbent + ever_switchedparty 
                             + I(Province)
                             , data = dnonrule[which(dnonrule$PR == 1), ], family = binomial(link = "logit"))
munmod.nonrule.obj.pr.c <- coeftest(munmod.nonrule.obj.pr, vcov = cluster.vcov(munmod.nonrule.obj.pr, dnonrule$cat_b[which(dnonrule$PR ==1)]))
# Create table
labs <- c("Service Rating Index", "Service Coverage $\\Delta$", 
          "Formal Housing $\\Delta$", "Log Population $\\Delta$",
          "Ethnicity - African (\\%)", "Post Secondary (\\%)", "Win Margin (2011)",
          "Years Incumbent", "Switched Party")
nmunis <- length(unique(drule$cat_b))
nmunis_anc <- length(unique(drule_ANC$cat_b))
stargazer(
          munmod.all.obj.ward.c,
          munmod.alla.obj.ward.c,
          munmod.nonrule.obj.ward.c,
          munmod.all.obj.pr.c,
          munmod.alla.obj.pr.c,
          munmod.nonrule.obj.pr.c,
          no.space = TRUE,
          covariate.labels = labs,
          column.labels = rep(c("All Ruling", "ANC", "Non-Ruling"), 2),
          omit = c("Province", "Constant"),
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = FALSE,
          add.lines = list(
            c("N Councillors", nobs(munmod.all.obj.ward), nobs(munmod.alla.obj.ward), nobs(munmod.nonrule.obj.ward),
            nobs(munmod.all.obj.pr), nobs(munmod.alla.obj.pr), nobs(munmod.nonrule.obj.pr)),
            c("N Municipalities", rep(c(nmunis, nmunis_anc, nmunis), 2)),
            c("Province FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")
            ),
          font.size = "footnotesize",
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes = "*$p < 0.05$",
          notes.label = "RSE clustered at municipal level",
          title = "Determinants of Councillor Renomination - All Municipalities (Logit), Ward vs. PR",
          label = "tab:munmain_wardvsprupdated"
          ,
          out = "munmain_wardvspr_updated.tex"
          )

## Appendix Table 13: Councillor Renomination and Political Competition - All Municipalities (Logit), Ward vs. PR ##

# Column 1: Service rating and win margin interaction - all ruling (WARD ONLY)
munmod.allx.ward <- glm(renom ~
                          service_rating_index2*win_margin
                        + service_pct_change_index1 
                        + formaldwell_pct_delta  
                        + logpopchange 
                        + black_avg 
                        + post_secondary_avg
                        + years_incumbent 
                        + ever_switchedparty 
                        + I(Province),
                        data = drule_ward, family = binomial(link = "logit"))
munmod.allx.ward.c <- coeftest(munmod.allx.ward, vcov = cluster.vcov(munmod.allx.ward, cluster = drule_ward$cat_b))
# Column 2: Service rating and win margin interaction - ANC only (WARD ONLY)
munmod.ancx.ward <- glm(renom ~ 
                          service_rating_index2*win_margin 
                        + service_pct_change_index1 
                        + formaldwell_pct_delta  
                        + logpopchange 
                        + black_avg 
                        + post_secondary_avg
                        + years_incumbent 
                        + ever_switchedparty 
                        + I(Province), 
                        data = drule_ANC[which(drule_ANC$Ward == 1), ], family = binomial(link = "logit"))
munmod.ancx.ward.c <- coeftest(munmod.ancx.ward, vcov = cluster.vcov(munmod.ancx.ward, cluster = drule_ANC$cat_b[which(drule_ANC$Ward == 1)]))
# Column 3: Service rating and win margin interaction - non-ruling (WARD ONLY)
munmod.nonrulex.ward <- glm(renom ~
                              service_rating_index2*win_margin 
                            + service_pct_change_index1  
                            + formaldwell_pct_delta  
                            + logpopchange 
                            + black_avg 
                            + post_secondary_avg
                            + years_incumbent 
                            + ever_switchedparty 
                            + I(Province),
                            data = dnonrule[which(dnonrule$Ward == 1), ], family = binomial(link = "logit"))
munmod.nonrulex.ward.c <- coeftest(munmod.nonrulex.ward, vcov = cluster.vcov(munmod.nonrulex.ward, cluster = dnonrule$cat_b[which(dnonrule$Ward == 1)]))
# Column 4: Service rating and win margin interaction - all ruling (PR ONLY)
munmod.allx.pr <- glm(renom ~ 
                        service_rating_index2*win_margin 
                      + service_pct_change_index1 
                      + formaldwell_pct_delta  
                      + logpopchange
                      + black_avg 
                      + post_secondary_avg
                      + years_incumbent 
                      + ever_switchedparty 
                      + I(Province)
                      , data = drule_pr, family = binomial(link = "logit"))
munmod.allx.pr.c <- coeftest(munmod.allx.pr, vcov = cluster.vcov(munmod.allx.pr, cluster = drule_pr$cat_b))
# Column 5: Service rating and win margin interaction - ANC only (PR ONLY)
munmod.ancx.pr <- glm(renom ~ 
                        service_rating_index2*win_margin 
                      + service_pct_change_index1 
                      + formaldwell_pct_delta  
                      + logpopchange
                      + black_avg 
                      + post_secondary_avg
                      + years_incumbent 
                      + ever_switchedparty 
                      + I(Province)
                      , data = drule_ANC[which(drule_ANC$PR == 1), ], family = binomial(link = "logit"))
munmod.ancx.pr.c <- coeftest(munmod.ancx.pr, vcov = cluster.vcov(munmod.ancx.pr, cluster = drule_ANC$cat_b[which(drule_ANC$PR == 1)]))
# Column 6: Service rating and win margin interaction - non-ruling (PR ONLY)
munmod.nonrulex.pr <- glm(renom ~ 
                            service_rating_index2*win_margin 
                          + service_pct_change_index1  
                          + formaldwell_pct_delta  
                          + logpopchange 
                          + black_avg 
                          + post_secondary_avg
                          + years_incumbent
                          + ever_switchedparty 
                          + I(Province)
                          , data = dnonrule[which(dnonrule$PR == 1), ], family = binomial(link = "logit"))
munmod.nonrulex.pr.c <- coeftest(munmod.nonrulex.pr, vcov = cluster.vcov(munmod.nonrulex.pr, cluster = dnonrule$cat_b[which(dnonrule$PR == 1)]))
# Create table
labs <- c("Svc. Rating Index x Win Margin", 
          "Service Rating Index", "Win Margin (2011)", "Service Coverage $\\Delta$",
          "Formal Housing $\\Delta$", "Log Population $\\Delta$",
          "Ethnicity - African (\\%)", "Post Secondary (\\%)",
          "Years Incumbent", "Switched Party")          
ncounc <- c(nobs(munmod.allx.ward), nobs(munmod.ancx.ward), nobs(munmod.nonrulex.ward),
            nobs(munmod.allx.pr), nobs(munmod.ancx.pr), nobs(munmod.nonrulex.pr))
nmunis <- length(unique(drule$cat_b))
nmunis_anc <- length(unique(drule_ANC$cat_b))
stargazer(munmod.allx.ward.c, munmod.ancx.ward.c, munmod.nonrulex.ward.c,
          munmod.allx.pr.c, munmod.ancx.pr.c, munmod.nonrulex.pr.c,
          no.space = TRUE,
          column.sep.width = "0.1pt",
          covariate.labels = labs,
          column.labels = rep(c("All Ruling", "ANC", "Non-Ruling"), 2),
          omit = c("Province", "Constant"),
          order = c(18, c(2:10)), 
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = FALSE,
          font.size = "footnotesize",
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes = "*$p<0.05$",
          add.lines = list(
            c("N Councillors", ncounc),
            c("N Municipalities", rep(c(nmunis, nmunis_anc, nmunis), 2)),
            c("Province FE", rep("Yes", 6))
          ),
          notes.label = "RSE clustered at municipal level",
          title = "Councillor Renomination and Political Competition - All Municipalities (Logit), Ward vs. PR",
          label = "tab:muncompetition_wardvsprupdated"
          ,
          out = "muncompetition_wardvspr_updated.tex"
          )

## Appendix Table 14: Councillor Renomination and Political Competition (Binary) in Gauteng (Logit) ##

# create binary measure of competitiveness - less than 35% vote margin
gcro_noBE$competitive_bin25 <- ifelse(gcro_noBE$win_margin2011 < 0.25, 1, 0)
gcro_noBE$competitive_bin35 <- ifelse(gcro_noBE$win_margin2011 < 0.35, 1, 0)
gcro_ANC <- gcro_noBE[gcro_noBE$Party.Name == "AFRICAN NATIONAL CONGRESS", ]
gcro_DA <- gcro_noBE[gcro_noBE$Party.Name == "DEMOCRATIC ALLIANCE", ]

# Column 1: All wards, binary vote margin
all.modxbin25 <- glm(renom ~ 
                     local_satisfied_15*competitive_bin25
                   + subj_sp_index15 
                   + obj_sp_delta_index_1315
                   + formal_housing_delta_1315 
                   + unemp_delta_1315_reverse
                   + high_income_avg 
                   + post_secondary_avg 
                   + black_avg
                   + ce_index15 
                   + logpop2011
                   + years_incumbent 
                   + ever_switchedparty 
                   + I(Municipality),
                   data = gcro_noBE, family = binomial(link = "logit"))
# Column 2: ANC wards, binary vote margin
anc.modxbin25 <- glm(renom ~ 
                     local_satisfied_15*competitive_bin25
                   + subj_sp_index15
                   + obj_sp_delta_index_1315
                   + formal_housing_delta_1315
                   + unemp_delta_1315_reverse
                   + high_income_avg
                   + post_secondary_avg
                   + black_avg
                   + ce_index15
                   + logpop2011
                   + years_incumbent
                   + ever_switchedparty 
                   + I(Municipality),
                   data = gcro_ANC, family = binomial(link = "logit"))
# Column 3: DA wards, binary vote margin
da.modxbin25 <- glm(renom ~ 
                    local_satisfied_15*competitive_bin25
                  + subj_sp_index15
                  + obj_sp_delta_index_1315
                  + formal_housing_delta_1315
                  + unemp_delta_1315_reverse
                  + subj_sp_index15 
                  + high_income_avg
                  + post_secondary_avg
                  + black_avg
                  + ce_index15
                  + logpop2011
                  + years_incumbent
                  + ever_switchedparty 
                  + I(Municipality),
                  data = gcro_DA, family = binomial(link = "logit"))
# Column 4: All wards, binary vote margin
all.modxbin35 <- glm(renom ~ 
                     local_satisfied_15*competitive_bin35
                   + subj_sp_index15
                   + obj_sp_delta_index_1315
                   + formal_housing_delta_1315
                   + unemp_delta_1315_reverse
                   + high_income_avg
                   + post_secondary_avg
                   + black_avg
                   + ce_index15
                   + logpop2011
                   + years_incumbent
                   + ever_switchedparty 
                   + I(Municipality),
                   data = gcro_noBE, family = binomial(link = "logit"))
# Column 5: ANC wards, binary vote margin
anc.modxbin35 <- glm(renom ~ 
                     local_satisfied_15*competitive_bin35
                   + subj_sp_index15
                   + obj_sp_delta_index_1315
                   + formal_housing_delta_1315
                   + unemp_delta_1315_reverse
                   + high_income_avg
                   + post_secondary_avg
                   + black_avg
                   + ce_index15
                   + logpop2011
                   + years_incumbent
                   + ever_switchedparty 
                   + I(Municipality),
                   data = gcro_ANC, family = binomial(link = "logit"))
# Column 6: DA wards, binary vote margin
da.modxbin35 <- glm(renom ~ 
                    local_satisfied_15*competitive_bin35
                  + subj_sp_index15
                  + obj_sp_delta_index_1315
                  + formal_housing_delta_1315
                  + unemp_delta_1315_reverse
                  + subj_sp_index15 
                  + high_income_avg
                  + post_secondary_avg
                  + black_avg
                  + ce_index15
                  + logpop2011
                  + years_incumbent
                  + ever_switchedparty 
                  + I(Municipality),
                  data = gcro_DA, family = binomial(link = "logit"))
# Create table
labs <- c("Councillor Satisfaction x Competitive (25\\%)",
          "Councillor Satisfaction x Competitive (35\\%)",
          "Councillor Satisfaction",
          "Competitive (25\\%)",
          "Competitive (35\\%)",
          "Service Satisfaction Index",
          "Service Coverage $\\Delta$", "Formal Housing $\\Delta$",
          "Employment $\\Delta$", "High Income (\\%)", "Post Secondary (\\%)",
          "Ethnicity - African (\\%)", "Civic Engagement",
          "Log Population (2011)", "Years Incumbent", "Switched Party")
stargazer(all.modxbin25,
          anc.modxbin25,
          da.modxbin25,
          all.modxbin35,
          anc.modxbin35,
          da.modxbin35,
          covariate.labels = labs,
          no.space = T,
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = F,
          font.size = "small",
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes.label = "",
          notes = "*$p<0.05$",
          title = "Councillor Renomination and Political Competition (Binary) in Gauteng (Logit)",
          omit = c("I\\(", "Constant"),
          omit.stat = c("ll", "aic"),
          order = c(24,25, 1:14),
          add.lines = list(c("Municipality FE", rep("Yes", 6))),
          column.labels = rep(c("All", "ANC", "DA"), 2),
          type = "latex",
          out = "wardcompetitionbin.tex",
          label = "tab:wardcompetitionbin")

## Appendix Table 15: Determinants of Councillor Renomination in Gauteng with Baseline Levels (Logit) ##

# Column 1: All Wards
all.mod.b <- glm(renom ~ 
                   local_satisfied_15
                 + subj_sp_index15
                 + obj_sp_delta_index_1315
                 + obj_sp_index13
                 + formal_housing_delta_1315 
                 + unemp_delta_1315_reverse
                 + high_income_avg 
                 + post_secondary_avg
                 + black_avg 
                 + ce_index15 
                 + logpop2011 
                 + win_margin2011
                 + years_incumbent 
                 + ever_switchedparty 
                 + I(Municipality),
                 data = gcro_noBE, family = binomial(link = "logit"))
# Column 2: ANC Only
anc.mod.b <- glm(renom ~ 
                   local_satisfied_15
                 + subj_sp_index15
                 + obj_sp_delta_index_1315
                 + obj_sp_index13
                 + formal_housing_delta_1315 
                 + unemp_delta_1315_reverse
                 + high_income_avg 
                 + post_secondary_avg
                 + black_avg 
                 + ce_index15 
                 + logpop2011 
                 + win_margin2011
                 + years_incumbent 
                 + ever_switchedparty 
                 + I(Municipality),
                 data = gcro_ANC, family = binomial(link = "logit"))
# Column 3: DA Only
da.mod.b <- glm(renom ~ 
                   local_satisfied_15
                 + subj_sp_index15
                 + obj_sp_delta_index_1315
                 + obj_sp_index13
                 + formal_housing_delta_1315 
                 + unemp_delta_1315_reverse
                 + high_income_avg 
                 + post_secondary_avg
                 + black_avg
                 + ce_index15
                 + logpop2011
                 + win_margin2011
                 + years_incumbent
                 + ever_switchedparty 
                 + I(Municipality),
                 data = gcro_DA, family = binomial(link = "logit"))
# Column 4: All Wards - interaction
all.mod.bx <- glm(renom ~ 
                   local_satisfied_15
                 + subj_sp_index15
                 + obj_sp_delta_index_1315*obj_sp_index13
                 + formal_housing_delta_1315 
                 + unemp_delta_1315_reverse
                 + high_income_avg 
                 + post_secondary_avg
                 + black_avg
                 + ce_index15
                 + logpop2011
                 + win_margin2011
                 + years_incumbent
                 + ever_switchedparty 
                 + I(Municipality),
                 data = gcro_noBE, family = binomial(link = "logit"))
# Column 5: ANC Only - interaction
anc.mod.bx <- glm(renom ~ 
                    local_satisfied_15
                  + subj_sp_index15
                  + obj_sp_delta_index_1315*obj_sp_index13
                  + formal_housing_delta_1315 
                  + unemp_delta_1315_reverse
                  + high_income_avg 
                  + post_secondary_avg
                  + black_avg
                  + ce_index15
                  + logpop2011
                  + win_margin2011
                  + years_incumbent
                  + ever_switchedparty 
                  + I(Municipality),
                  data = gcro_ANC, family = binomial(link = "logit"))
# Column 6: DA Only - interaction
da.mod.bx <- glm(renom ~ 
                    local_satisfied_15
                  + subj_sp_index15
                  + obj_sp_delta_index_1315*obj_sp_index13
                  + formal_housing_delta_1315 
                  + unemp_delta_1315_reverse
                  + high_income_avg 
                  + post_secondary_avg
                  + black_avg
                  + ce_index15
                  + logpop2011
                  + win_margin2011
                  + years_incumbent
                  + ever_switchedparty 
                  + I(Municipality),
                  data = gcro_DA, family = binomial(link = "logit"))
# Create table
labs <- c("Councillor Satisfaction",
          "Service Satisfaction Index",
          "Service Coverage $\\Delta$", "Service Coverage (2013)",
          "Formal Housing $\\Delta$",
          "Employment $\\Delta$", "High Income (\\%)", "Post Secondary (\\%)",
          "Ethnicity - African (\\%)", "Civic Engagement",
          "Log Population (2011)", "Win Margin (2011)", "Years Incumbent", "Switched Party",
          "Service $\\Delta$ X Service Coverage (2013)")
stargazer(all.mod.b, anc.mod.b, da.mod.b,
          all.mod.bx, anc.mod.bx, da.mod.bx,
          no.space = TRUE,
          covariate.labels = labs,
          column.labels = rep(c("All", "ANC", "DA"), 2),
          omit = c(" - ", "Constant"),
          omit.stat = c("ll", "aic"),
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = FALSE,
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes = "*$p<0.05$",
          font.size = "footnotesize",
          title = "Determinants of Councillor Renomination in Gauteng with Baseline Levels (Logit)",
          add.lines = list(c("Municipality FE", rep("Yes", 6))),
          label = "tab:wardmainbline"
         ,
         out = "wardmainbline.tex"
         )

## Appendix Table 16: Determinants of Councillor Renomination in Gauteng - Excluding Fully Serviced Wards (Logistic Regression) ##

# Subset this to ANC only
gcro_nfs_ANC <- gcro_nfs[which(gcro_nfs$Party.Name == "AFRICAN NATIONAL CONGRESS"), ]

# Column 1: All Wards
all.mod.nfs <- glm(renom ~ 
                     local_satisfied_15
                   + subj_sp_index15
                   + obj_sp_delta_index_1315 
                   + formal_housing_delta_1315
                   + unemp_delta_1315_reverse
                   + high_income_avg
                   + post_secondary_avg
                   + black_avg
                   + ce_index15
                   + logpop2011
                   + win_margin2011
                   + years_incumbent
                   + ever_switchedparty 
                   + I(Municipality),
                   data = gcro_nfs, family = binomial(link = "logit"))
# Column 2: ANC Only
anc.mod.nfs <- glm(renom ~ 
                     local_satisfied_15
                   + subj_sp_index15
                   + obj_sp_delta_index_1315 
                   + formal_housing_delta_1315
                   + unemp_delta_1315_reverse
                   + high_income_avg
                   + post_secondary_avg
                   + black_avg
                   + ce_index15
                   + logpop2011
                   + win_margin2011
                   + years_incumbent
                   + ever_switchedparty 
                   + I(Municipality),
                   data = gcro_nfs_ANC, family = binomial(link = "logit"))
labs <- c( "Councillor Satisfaction",
          "Service Satisfaction Index",
          "Service Coverage $\\Delta$",
          "Formal Housing $\\Delta$",
          "Employment $\\Delta$", "High Income (\\%)", "Post Secondary (\\%)",
          "Ethnicity - African (\\%)", "Civic Engagement",
          "Log Population (2011)", 
          "Win Margin (2011)",
          "Years Incumbent", "Switched Party")
stargazer(all.mod.nfs,
          anc.mod.nfs, 
          no.space = TRUE,
          covariate.labels = labs,
          column.labels = c("All", "ANC"),
          omit = c(" - ", "Constant"),
          omit.stat = c("ll", "aic"),
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = FALSE,
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes.label = "",
          notes = "*$p<0.05$",
          title = "Determinants of Councillor Renomination in Gauteng - Excluding Fully Serviced Wards (Logistic Regression)",
          add.lines = list(c("Municipality FE", "Yes", "Yes")),
          label = "tab:wardmainnfs"
          ,
          out = "wardmainnfs.tex"
          )

## Appendix Table 17: Determinants of Councillor Renomination with Baseline Levels - All Municipalities (Logit) ##

# Column 1: All ruling
munmod.all.b <- glm(renom ~ 
                      service_rating_index2
                    + service_pct_change_index1 
                    + service_2011_index1
                    + formaldwell_pct_delta  
                    + logpopchange
                    + black_avg
                    + post_secondary_avg
                    + win_margin
                    + Ward
                    + years_incumbent
                    + ever_switchedparty
                    + I(Province), data = drule_top10, family = binomial(link = "logit"))
munmod.all.b.c <- coeftest(munmod.all.b, vcov = cluster.vcov(munmod.all.b, drule_top10$cat_b))
# Column 2: ANC only
munmod.anc.b <- glm(renom ~ 
                      service_rating_index2
                    + service_pct_change_index1 
                    + service_2011_index1
                    + formaldwell_pct_delta  
                    + logpopchange
                    + black_avg
                    + post_secondary_avg
                    + win_margin
                    + Ward
                    + years_incumbent
                    + ever_switchedparty
                    + I(Province), data = drule_ANC_top10, family = binomial(link = "logit"))
munmod.anc.b.c <- coeftest(munmod.anc.b, vcov = cluster.vcov(munmod.anc.b, drule_ANC_top10$cat_b))
# Column 3: Non-ruling only
munmod.nonrule.b <- glm(renom ~ 
                          service_rating_index2
                        + service_pct_change_index1 
                        + service_2011_index1
                        + formaldwell_pct_delta  
                        + logpopchange
                        + black_avg
                        + post_secondary_avg
                        + win_margin
                        + Ward
                        + years_incumbent
                        + ever_switchedparty
                        + I(Province), data = dnonrule_top10, family = binomial(link = "logit"))
munmod.nonrule.b.c <- coeftest(munmod.nonrule.b, vcov = cluster.vcov(munmod.nonrule.b, dnonrule_top10$cat_b))
# Column 4: All ruling - interaction
munmod.all.bx <- glm(renom ~ 
                       service_rating_index2 
                     + service_pct_change_index1*service_2011_index1
                     + formaldwell_pct_delta  
                     + logpopchange
                     + black_avg
                     + post_secondary_avg
                     + win_margin
                     + Ward
                     + years_incumbent
                     + ever_switchedparty
                     + I(Province),
                     data = drule_top10, family = binomial(link = "logit"))
munmod.all.bx.c <- coeftest(munmod.all.bx, vcov = cluster.vcov(munmod.all.bx, drule_top10$cat_b))
# Column 5: ANC only - interaction
munmod.anc.bx <- glm(renom ~ 
                       service_rating_index2 
                     + service_pct_change_index1*service_2011_index1
                     + formaldwell_pct_delta  
                     + logpopchange
                     + black_avg
                     + post_secondary_avg
                     + win_margin
                     + Ward
                     + years_incumbent
                     + ever_switchedparty
                     + I(Province),
                     data = drule_ANC_top10, family = binomial(link = "logit"))
munmod.anc.bx.c <- coeftest(munmod.anc.bx, vcov = cluster.vcov(munmod.anc.bx, drule_ANC_top10$cat_b))
# Column 6: Non-ruling only - interaction
munmod.nonrule.bx <- glm(renom ~ 
                           service_rating_index2 
                         + service_pct_change_index1*service_2011_index1
                         + formaldwell_pct_delta  
                         + logpopchange
                         + black_avg
                         + post_secondary_avg
                         + win_margin
                         + Ward
                         + years_incumbent
                         + ever_switchedparty
                         + I(Province),
                         data = dnonrule_top10, family = binomial(link = "logit"))
munmod.nonrule.bx.c <- coeftest(munmod.nonrule.bx, vcov = cluster.vcov(munmod.nonrule.bx, dnonrule_top10$cat_b))
# Create table
labs <- c( "Service Rating Index",
           "Service Coverage $\\Delta$",
           "Service Coverage Index (2011)",
           "Formal Housing $\\Delta$", "Log Population $\\Delta$",
           "Ethnicity - African (\\%)", "Post Secondary (\\%)",
           "Win Margin (2011)",
           "Years Incumbent", "Switched Party", "Service $\\Delta$ X Service Coverage (2011)")
nmunis <- length(unique(drule_top10$cat_b))
nmunis_anc <- length(unique(drule_ANC_top10$cat_b))
stargazer(munmod.all.b.c, munmod.anc.b.c, munmod.nonrule.b.c,
          munmod.all.bx.c, munmod.anc.bx.c, munmod.nonrule.bx.c,
          no.space = TRUE,
          covariate.labels = labs,
          font.size = "small",
          dep.var.caption = "Dependent variable: Renomination",
          dep.var.labels.include = FALSE,
          star.char = c("*"),
          star.cutoffs = c(0.05),
          notes.append = F,
          notes = "*$p<0.05$",
          column.labels = rep(c("All Ruling", "ANC", "Non-Ruling"), 2),
          omit = c("Province", "Constant"),
          add.lines = list(
            c("N Councillors", rep(c(nobs(munmod.all.b), nobs(munmod.anc.b), nobs(munmod.nonrule.b)), 2)),
            c("N Municipalities", rep(c(nmunis, nmunis_anc, nmunis), 2)),
            c("Province FE", rep("Yes", 6))
          ),
          notes.label = "RSE clustered at municipal level",
          title = "Determinants of Councillor Renomination with Baseline Levels - Top 10 PR Councillors, All Municipalities (Logit)",
          label = "tab:munmainbline_top10"
          ,
          out = "munmainbline_top10.tex"
          )

## Appendix Table 18: Sensitivity Statistics (Gauteng) ##

# Column 3: All wards - councillor satisfaction, service satisfaction, objective services
all.mod.obj.ols <- lm(renom ~ 
                        local_satisfied_15 
                      + subj_sp_index15
                      + obj_sp_delta_index_1315
                      + formal_housing_delta_1315
                      + unemp_delta_1315_reverse
                      + high_income_avg
                      + post_secondary_avg
                      + black_avg
                      + ce_index15
                      + logpop2011
                      + win_margin2011
                      + years_incumbent
                      + ever_switchedparty 
                      + I(Municipality),
                      data = gcro_noBE)

sens.all <- sensemakr(all.mod.obj.ols, "local_satisfied_15", c("subj_sp_index15", 
                                                               "obj_sp_delta_index_1315",
                                                               "formal_housing_delta_1315",
                                                               "unemp_delta_1315_reverse",
                                                               "high_income_avg",
                                                               "post_secondary_avg",
                                                               "black_avg",
                                                               "ce_index15",
                                                               "logpop2011",
                                                               "win_margin2011",
                                                               "years_incumbent",
                                                               "ever_switchedparty"
), kd = c(1, 2, 3))

# Create table
sens.table.print <- print(ovb_minimal_reporting(sens.all))
save(
     file = "sens_table.tex",
     sens.table.print)

## Appendix Table 19: Omitted Variable Bias Bounds (Sensitivity) ##
sense.covs.x <- xtable(sens.all$bounds[, -c(4, 8, 9)], label = "tab:senscovs",
                       caption = "Omitted Variable Bias Bounds (Sensitivity)")
print(sense.covs.x,
      file = "sens_covs.tex",
      caption.placement = "top",
      include.rownames = F)

## Appendix Table 20: Objective Service Provision and Councillor Satisfaction - Afrobarometer R6 (Individual) ##

### PLEASE CONTACT AUTHORS FOR INDIVIDUAL-LEVEL DATA REQUIRED TO PRODUCE THIS TABLE ###

# # Column 1: Municipality Service Change
# ab_mod1 <- lm(localcouncilorperformance_subj ~
#                 mun_service_pct_change_index1
#               + female 
#               + race_black 
#               + education 
#               + havejob 
#               + owngoods.index 
#               + provFE,
#               data = ab)
# ab_mod1SE <- coeftest(ab_mod1, vcov = cluster.vcov(ab_mod1, ab$cat_b))
# # Column 2: Municipality Service Level
# ab_mod2 <- lm(localcouncilorperformance_subj ~
#                 mun_service_2016_index1
#               + female 
#               + race_black 
#               + education 
#               + havejob 
#               + owngoods.index 
#               + provFE,
#               data = ab)
# ab_mod2SE <- coeftest(ab_mod2, vcov = cluster.vcov(ab_mod2, ab$cat_b))
# # Column 3: Enumeration Area Service Level (Index 1: electricity, pipedwater, sewage, pavedroad)
# ab_mod3 <- lm(localcouncilorperformance_subj ~
#                 ea_serviceindex1
#               + female 
#               + race_black 
#               + education 
#               + havejob 
#               + owngoods.index 
#               + provFE,
#               data = ab)
# ab_mod3SE <- coeftest(ab_mod3, vcov = cluster.vcov(ab_mod3, ab$eanumb))
# # Column 4: Enumeration Area Service Level (Index 2: Post office, school, police station, healthclinic, cell service)
# ab_mod4 <- lm(localcouncilorperformance_subj ~
#                 ea_serviceindex2
#               + female 
#               + race_black 
#               + education 
#               + havejob 
#               + owngoods.index 
#               + provFE,
#               data = ab)
# ab_mod4SE <- coeftest(ab_mod4, vcov = cluster.vcov(ab_mod4, ab$eanumb))
# # Column 5: Household Service level (Index of water, toilet, electricity)
# ab$obj_sp_index_indiv <- ab$waterinsidehouse + ab$toiletinsidehouse + ab$electricityavailable_dummy
# ab_mod5 <- lm(localcouncilorperformance_subj ~
#                 obj_sp_index_indiv
#               + female
#               + race_black
#               + education
#               + havejob
#               + owngoods.index
#               + provFE,
#               data = ab)
# ab_mod5SE <- coeftest(ab_mod5, vcov = cluster.vcov(ab_mod5, ab$eanumb))
# # Column 6: Additional covariates
# ab_mod6 <- lm(localcouncilorperformance_subj ~
#                 obj_sp_index_indiv
#               + livingconditions_subj
#               + thinklocalcouncilorslisten
#               + contacted_councillor
#               + viewascorrupt_localcouncilors
#               + trust_rulingparty
#               + female
#               + race_black
#               + education
#               + havejob
#               + owngoods.index
#               + provFE,
#               data = ab)
# ab_mod6SE <- coeftest(ab_mod6, vcov = cluster.vcov(ab_mod6, ab$cat_b))
# # Create table
# labs <- c("Municipal Services $\\Delta$", "Municipal Services 2016", 
#           "Enumeration Area Services: Index 1", "Enumeration Area Services: Index 2",
#           "Household Service Index", "Evaluation of Living Conditions",
#           "Think Councillors Listen", "Contacted Councillor", "Think Councillors Corrupt",
#           "Trust Ruling Party",
#           "Female", "Ethnicity - African", "Education", "Employed",
#           "Asset Index")
# numobs <- c(nobs(ab_mod1), nobs(ab_mod2), nobs(ab_mod3), nobs(ab_mod4),
#             nobs(ab_mod5), nobs(ab_mod6))
# stargazer(ab_mod1SE, ab_mod2SE, ab_mod3SE, ab_mod4SE, ab_mod5, 
#           ab_mod6SE,
#           covariate.labels = labs,
#           no.space = TRUE,
#           column.sep.width = "1pt",
#           omit = c("provFE", "Constant"),
#           omit.stat = c("rsq", "adj.rsq", "ser", "f", "n"),
#           model.names = FALSE,
#           dep.var.caption = "DV: Approval of Councillor Performance",
#           dep.var.labels.include = FALSE, 
#           add.lines = list(
#             c("Observations", numobs),
#             c("Province FE", rep("Yes", 7)),
#             c("Cluster", c("Mun", "Mun", "EA", "EA", "", "Mun"))),
#           star.char = c("*"),
#           star.cutoffs = c(0.05),
#           notes.append = F,
#           notes.label = "",
#           notes = "*$p<0.05$",
#           title = "Objective Service Provision and Councilor Satisfaction - Afrobarometer R6 (Individual)",
#           label = "tab:abindiv"
#           ,
#           out = "abindiv.tex")

## Appendix Table 21: Objective Service Provision and Local Councillor Satisfaction in Gauteng (Individual) ##

### PLEASE CONTACT AUTHORS FOR INDIVIDUAL-LEVEL DATA REQUIRED TO PRODUCE THIS TABLE ###

# mod.obj.change <- lm(satisfied_lc ~ 
#                        obj_sp_delta_index_1315
#                      + refuse_removal_weekly
#                      + female 
#                      + edu_matric
#                      + race_african
#                      + income
#                      + formal_dwelling
#                      + unemp, data = gind)
# summary(mod.obj.change)
# mod.obj.change.c <- coeftest(mod.obj.change, vcov = cluster.vcov(mod.obj.change, gind$WardNumber))
# mod.obj.change.c
# # Column 2: Piped water
# mod.obj.water <- lm(satisfied_lc ~ 
#                       piped_water_home
#                     + female 
#                     + edu_matric
#                     + race_african
#                     + income 
#                     + formal_dwelling
#                     + unemp, data = gind)
# # Columnn 3: Toilet
# mod.obj.toilet <- lm(satisfied_lc ~ 
#                        flush_toilet
#                      + female 
#                      + edu_matric
#                      + race_african
#                      + income
#                      + formal_dwelling
#                      + unemp, data = gind)
# # Column 4: Electricity
# mod.obj.elec <- lm(satisfied_lc ~ 
#                      electricity_light
#                    + female 
#                    + edu_matric
#                    + race_african
#                    + income 
#                    + formal_dwelling
#                    + unemp, data = gind)
# # Column 5: Refuse removal
# mod.obj.refuse <- lm(satisfied_lc ~ 
#                        refuse_removal_weekly
#                      + female 
#                      + edu_matric
#                      + race_african
#                      + income 
#                      + formal_dwelling
#                      + unemp, data = gind)
# # Column 6: Objective services index
# mod.obj.index <- lm(satisfied_lc ~ 
#                       obj_sp_index_indiv
#                     + female 
#                     + edu_matric
#                     + race_african
#                     + income 
#                     + formal_dwelling
#                     + unemp, data = gind)
# 
# # Create table
# labs <- c("Service Coverage $\\Delta$", "Piped Water", "Flush Toilet", "Electricity", "Refuse Removal", "Service Index (2015)",
#           "Female", "Education (Matric)", "Ethnicity - African", "Income", "Formal Dwelling",
#           "Unemployed")
# numobs <- c(nobs(mod.obj.change), nobs(mod.obj.water), nobs(mod.obj.toilet), nobs(mod.obj.elec), nobs(mod.obj.refuse), nobs(mod.obj.index))
# stargazer(mod.obj.change.c, mod.obj.water, mod.obj.toilet, mod.obj.elec, mod.obj.refuse,
#           mod.obj.index,
#           model.names = F,
#           covariate.labels = labs,
#           no.space = TRUE,
#           omit = c("Constant"),
#           omit.stat = c("f", "ser", "n", "rsq", "adj.rsq"),
#           star.char = c("*"),
#           star.cutoffs = c(0.05),
#           notes.append = F,
#           notes.label = "",
#           notes = "*$p<0.05$",
#           dep.var.caption = "DV: Satisfaction with Local Councillor Performance",
#           dep.var.labels.include = FALSE,
#           add.lines = list(
#             c("Observations", numobs),
#             c("Cluster", "Ward", rep("None", 5))),
#           title = "Objective Service Provision and Local Councillor Satisfaction in Gauteng (Individual)",
#           label = "tab:gindiv",
#           out = "gindiv.tex")
# 


#### APPENDIX FIGURES ####

## Appendix Figure 1: Distribution of 2016 promotion outcomes for PR councillors elected in 2011, all municipalities in South Africa. ##

pdf("promotion_hist_muni_listmove.pdf")
xx <- barplot(table(d_pr$pr_listorder_partypercentiledelta2num),
              names.arg = c("not nominated", "renom \n lower quartile",
                            "ward nomination",
                            "renom \n same quartile",
                            "renom \n higher quartile",
                            "nominated to \n higher office"),
              ylim = c(0, 2250),
              main = "2016 Nomination Outcomes \n 2011 PR Councillors (n = 2,287)",
              cex.main = 0.9,
              cex.names = 0.7)
text(x = xx,
     y = as.numeric(table(d_pr$pr_listorder_partypercentiledelta2num)),
     labels = as.numeric(table(d_pr$pr_listorder_partypercentiledelta2num)),
     pos = 3,
     cex = 0.8)
dev.off()

## Appendix Figure 2: Distribution of 2016 nomination outcomes for ward councillors elected in 2011 in Gauteng province. ##

pdf("promotion_hist_ward.pdf")
xx <- barplot(table(gcro_noBE$renom_type_num),
              names.arg = c("no nomination", "ward \n nomination",
                            "PR \n nomination", "nominated to \n higher office"),
              ylim = c(0, 240),
             # main = "2016 Nomination Outcomes \n Gauteng 2011 Ward Councillors",
              cex.main = 0.9,
              cex.names = 0.7)
text(x = xx,
     y = as.numeric(table(gcro_noBE$renom_type_num)),
     labels = as.numeric(table(gcro_noBE$renom_type_num)),
     pos = 3,
     cex = 0.8)
dev.off()

## Appendix Figure 3: Distribution of estimates of coefficients for Service Coverage Change and Councillor Satisfaction with clustered bootstrap. ##

### PLEASE CONTACT AUTHORS FOR INDIVIDUAL-LEVEL DATA REQUIRED FOR THIS FIGURE ###

# Load gcro survey data and supplementary datasets needed

# R = 200 #number of bootstrap simulations to run
# B = 22 #number of coefficients estimated by model
# results = matrix(nrow = B, ncol= R, data = NA) #storage matrix for bootstrap estimates
# 
# 
# ##Bootstrap for-loop : creates replicate datasets by sampling within wards with replacement, then stores model coefficients
# 
# for(b in 1:R){
#   
#   #2013
#   
#   boot.data = gcrosub13
#   boot.numbers <- NA
#   for(i in 1:length(wards13)){
#     #select ward
#     wardsubset <- gcrosub13[which(gcrosub13$Ward==wards13[i]),]
#     #sample respondents with replacement
#     samplenumbers <- sample(wardsubset$X, size=length(wardsubset$X), replace=TRUE)
#     boot.numbers = c(boot.numbers, samplenumbers)
#     #print(i)
#   }
#   boot.numbers <- boot.numbers[-1]
#   boot.numbers <- as.matrix(boot.numbers)
#   colnames(boot.numbers) <- "X"
#   #now create replicate dataset
#   replicate13 = merge(boot.numbers, boot.data, by="X")
#   
#   # Create ward-level variables with replicate dataset
#   
#   myvars <- c("water_source", "elec", "flush_toilet", "waste_disposal", "crime_victim_no",
#               "unemp", "formal_housing", "trad_housing",
#               "water_satisfied", "sanitation_satisfied", "waste_satisfied", "energy_satisfied",
#               "roads_satisfied", "safety_satisfied", "batho_pele", "national_satisfied",
#               "provincial_satisfied", "muni_satisfied", "high_income", "high_status",
#               "business_owner", "business_owner_formal", "post_secondary", "black",
#               "reg_voter", "plan_vote", "attend_meeting", "heard_idp", "participate_idp", "participate_club",
#               "participate_protest")
#   
#   replicate13_wardvars13 <- matrix(NA, nrow = length(wards13), ncol = length(myvars))
#   for(i in 1:length(myvars)){
#     myvar <- replicate13[, myvars[i]]
#     replicate13_wardvars13[, i] <- by(myvar, replicate13$Ward, mean, na.rm = T)
#   }
#   replicate13_wardvars13 <- data.frame(replicate13_wardvars13)
#   colnames(replicate13_wardvars13) <- paste(myvars, "_13", sep = "")
#   rownames(replicate13_wardvars13) <- wards13
#   head(replicate13_wardvars13)
#   
#   # 2015
#   
#   boot.data=gcrosub15
#   boot.numbers <- NA
#   for(i in 1:length(wards15)){
#     #select ward
#     wardsubset <- gcrosub15[which(gcrosub15$Ward==wards15[i]),]
#     #sample respondents with replacement
#     samplenumbers <- sample(wardsubset$X, size=length(wardsubset$X), replace=TRUE)
#     boot.numbers = c(boot.numbers, samplenumbers)
#     #print(i)
#   }
#   boot.numbers <- boot.numbers[-1]
#   boot.numbers <- as.matrix(boot.numbers)
#   colnames(boot.numbers) <- "X"
#   #now create replicate dataset
#   replicate15 = merge(boot.numbers, boot.data, by="X")
#   
#   myvars15 <- c(myvars, "local_satisfied", "know_lc")
#   
#   replicate15_wardvars15 <- matrix(NA, nrow = length(wards15), ncol = length(myvars15))
#   for(i in 1:length(myvars15)){
#     myvar <- replicate15[, myvars15[i]]
#     replicate15_wardvars15[, i] <- by(myvar, replicate15$WardNumber, mean, na.rm = T)
#   }
#   replicate15_wardvars15 <- data.frame(replicate15_wardvars15)
#   colnames(replicate15_wardvars15) <- paste(myvars15, "_15", sep = "")
#   rownames(replicate15_wardvars15) <- wards15
#   head(replicate15_wardvars15)
#   
#   all.equal(row.names(replicate13_wardvars13), row.names(replicate15_wardvars15))
#   
#   # combine
#   replicate_gcro_ward <- cbind(replicate13_wardvars13, replicate15_wardvars15)
#   
#   # Create indices and Standardize
#   zFun <- function(data, varnames){
#     z_vars <- matrix(NA, nrow = nrow(data), ncol = length(varnames))
#     for(i in 1:length(varnames)){
#       mycol <- data[, varnames[i]]
#       mycol_mean <- mean(mycol, na.rm = T)
#       mycol_sd <- sd(mycol, na.rm = T)
#       mycol_std <- (mycol - mycol_mean)/mycol_sd
#       z_vars[, i] <- mycol_std
#     }
#     z_vars <- data.frame(z_vars)
#     colnames(z_vars) <- paste(varnames, "_std", sep = "")
#     rownames(z_vars) <- rownames(data)
#     return(z_vars)
#   }
#   
#   replicate_gcro_ward_std <- zFun(replicate_gcro_ward, colnames(replicate_gcro_ward))
#   
#   # 2013
#   
#   # Piped water, flush toilet, weekly waste removal, and electricity
#   replicate_gcro_obj_sp_vars13 <- replicate_gcro_ward_std[, c("water_source_13_std", "flush_toilet_13_std", "waste_disposal_13_std", "elec_13_std")]
#   replicate_gcro_ward$obj_sp_index13 <- apply(replicate_gcro_obj_sp_vars13, 1, mean)
#   
#   # Piped water, flush toilet, weekly waste removal, and not crime victimization (original analysis)
#   replicate_gcro_obj_sp_vars13a <- replicate_gcro_ward_std[, c("water_source_13_std", "flush_toilet_13_std", "waste_disposal_13_std", "crime_victim_no_13_std")]
#   replicate_gcro_ward$obj_sp_index13a <- apply(replicate_gcro_obj_sp_vars13a, 1, mean)
#   
#   # 2015
#   
#   # Piped water, flush toilet, weekly waste removal, and electricity
#   replicate_gcro_obj_sp_vars15 <- replicate_gcro_ward_std[, c("water_source_15_std", "flush_toilet_15_std", "waste_disposal_15_std", "elec_15_std")]
#   replicate_gcro_ward$obj_sp_index15 <- apply(replicate_gcro_obj_sp_vars15, 1, mean)
#   
#   # Piped water, flush toilet, weekly waste removal, and not crime victimization (original analysis)
#   replicate_gcro_obj_sp_vars15a <- replicate_gcro_ward_std[, c("water_source_15_std", "flush_toilet_15_std", "waste_disposal_15_std", "crime_victim_no_15_std")]
#   replicate_gcro_ward$obj_sp_index15a <- apply(replicate_gcro_obj_sp_vars15a, 1, mean)
#   
#   # Subjective Service Provision Index
#   
#   # 2013
#   
#   # Only services we have objective measures for (water, sanitation, waste, energy)
#   replicate_gcro_subj_sp_vars13 <- replicate_gcro_ward_std[, c("water_satisfied_13_std", "sanitation_satisfied_13_std",
#                                                                "waste_satisfied_13_std", "energy_satisfied_13_std")]
#   replicate_gcro_ward$subj_sp_index13 <- apply(replicate_gcro_subj_sp_vars13, 1, mean)
#   
#   # Original (all service satisfaction measures)
#   replicate_gcro_subj_sp_vars13a <- replicate_gcro_ward_std[, c("water_satisfied_13_std", "sanitation_satisfied_13_std",
#                                                                 "waste_satisfied_13_std", "energy_satisfied_13_std", 
#                                                                 "roads_satisfied_13_std", "safety_satisfied_13_std")]
#   replicate_gcro_ward$subj_sp_index13a <- apply(replicate_gcro_subj_sp_vars13a, 1, mean)
#   
#   # 2015
#   
#   # Only services we have objective measures for (water, sanitation, waste, energy)
#   replicate_gcro_subj_sp_vars15 <- replicate_gcro_ward_std[, c("water_satisfied_15_std", "sanitation_satisfied_15_std",
#                                                                "waste_satisfied_15_std", "energy_satisfied_15_std")]
#   replicate_gcro_ward$subj_sp_index15 <- apply(replicate_gcro_subj_sp_vars15, 1, mean)
#   
#   # Original (all service satisfaction measures)
#   replicate_gcro_subj_sp_vars15a <- replicate_gcro_ward_std[, c("water_satisfied_15_std", "sanitation_satisfied_15_std",
#                                                                 "waste_satisfied_15_std", "energy_satisfied_15_std", 
#                                                                 "roads_satisfied_15_std", "safety_satisfied_15_std")]
#   replicate_gcro_ward$subj_sp_index15a <- apply(replicate_gcro_subj_sp_vars15a, 1, mean)
#   
#   # Civic Engagement Index
#   
#   replicate_gcro_ce_13 <- replicate_gcro_ward_std[, c("reg_voter_13_std", "plan_vote_13_std", "attend_meeting_13_std",
#                                                       "heard_idp_13_std", "participate_idp_13_std", "participate_club_13_std",
#                                                       "participate_protest_13_std")]
#   replicate_gcro_ward$ce_index13 <- apply(replicate_gcro_ce_13, 1, mean)
#   
#   replicate_gcro_ce_15 <- replicate_gcro_ward_std[, c("reg_voter_15_std", "plan_vote_15_std", "attend_meeting_15_std",
#                                                       "heard_idp_15_std", "participate_idp_15_std", "participate_club_15_std",
#                                                       "participate_protest_15_std")]
#   replicate_gcro_ward$ce_index15 <- apply(replicate_gcro_ce_15, 1, mean)
#   
#   # Black (all 3 years)
#   replicate_gcro_ward$black_avg <- apply(replicate_gcro_ward[, c("black_13", "black_15")], 1, mean, na.rm = T)
#   # Average civic education index (all 3 years)
#   replicate_gcro_ward$ce_avg <- apply(replicate_gcro_ward[, c("ce_index13", "ce_index15")], 1, mean, na.rm = T)
#   # Average pct high income (all 3 years)
#   replicate_gcro_ward$high_income_avg <- apply(replicate_gcro_ward[, c("high_income_13", "high_income_15")], 1, mean, na.rm = T)
#   # Average pct with post-secondary education
#   replicate_gcro_ward$post_secondary_avg <- apply(replicate_gcro_ward[, c("post_secondary_13", "post_secondary_15")], 1, mean, na.rm = T)
#   
#   # NEW INDICES STANDARDIZED BY MEAN CHANGE (DELTA INDEX)
#   
#   # Objective service provision
#   
#   # Change between 2013 and 2015
#   replicate_gcro_ward$water_source_delta_1315 <- replicate_gcro_ward$water_source_15 - replicate_gcro_ward$water_source_13
#   replicate_gcro_ward$flush_toilet_delta_1315 <- replicate_gcro_ward$flush_toilet_15 - replicate_gcro_ward$flush_toilet_13 
#   replicate_gcro_ward$waste_disposal_delta_1315 <- replicate_gcro_ward$waste_disposal_15 - replicate_gcro_ward$waste_disposal_13
#   replicate_gcro_ward$elec_delta_1315 <- replicate_gcro_ward$elec_15 - replicate_gcro_ward$elec_13
#   
#   # Standardizing the differences
#   replicate_obj_sp_deltas_std <- zFun(replicate_gcro_ward, c("water_source_delta_1315", "flush_toilet_delta_1315",
#                                                              "waste_disposal_delta_1315", "elec_delta_1315"
#   ))
#   # Mean of standardized differences
#   replicate_gcro_ward$obj_sp_delta_index_1315 <- apply(replicate_obj_sp_deltas_std[, c(1:4)], 1, mean)
#   
#   # Change in unemployment
#   replicate_gcro_ward$unemp_delta_1315 <- replicate_gcro_ward$unemp_15 - replicate_gcro_ward$unemp_13
#   
#   # Reverse (so positive change represents decrease in unemployment)
#   replicate_gcro_ward$unemp_delta_1315_reverse <- -1 * replicate_gcro_ward$unemp_delta_1315
#   
#   # Change in formal housing
#   replicate_gcro_ward$formal_housing_delta_1315 <- replicate_gcro_ward$formal_housing_15 - replicate_gcro_ward$formal_housing_13
#   
#   # Merge in municipality variable
#   replicate_gcro_ward$ward <- rownames(replicate_gcro_ward)
#   munmerg <- as.data.frame(as.matrix(cbind(as.character(g15$Municipality), g15$WardNumber)))
#   colnames(munmerg) <- c("muni", "ward")
#   replicate_gcro_ward <- merge(replicate_gcro_ward, munmerg, by="ward")
#   #replicate_gcro_ward$muni <- NA
#   #for(i in 1:nrow(replicate_gcro_ward)){
#   #  replicate_gcro_ward$muni[i] <- unique(as.character(g15$Municipality[which(g15$WardNumber == replicate_gcro_ward$ward[i])]))
#   #print(i)
#   #  }
#   
#   # Merge in populations data from 2011 census
#   cen$ward <- cen$ward_id
#   censub <- cen[, c("ward", "pop2011")]
#   censub_gt <- censub[which(censub$ward %in% replicate_gcro_ward$ward), ]
#   replicate_gcro_ward <- merge(replicate_gcro_ward, censub_gt, "ward")
#   
#   # Log population
#   replicate_gcro_ward$logpop2011 <- log(replicate_gcro_ward$pop2011)
#   
#   # Merge in By-election variables and additional variables
#   gcro_councillors <- read.csv("./Scratchpad/GCRO_analysis/gcro_councillors.csv")
#   add_vars <- as.data.frame(cbind(gcro_councillors$ward, gcro_councillors$BE, as.character(gcro_councillors$promotionoutcome2), gcro_councillors$ever_switchedmuni, gcro_councillors$win_margin2011))
#   colnames(add_vars) <- c("ward", "BE", "promotionoutcome2", "ever_switchedmuni", "win_margin2011")
#   replicate_gcro_ward <- merge(replicate_gcro_ward, add_vars, by = "ward")
#   
#   # Merge all GCRO variables with 2011 councillor list
#   replicate_merged <- merge(wcs_gt, replicate_gcro_ward, by = "ward")
#   
#   # Create subsets for analysis
#   replicate_gcro_noBE <- replicate_merged[which(replicate_merged$BE == 0),]
#   replicate_gcro_ANC <- replicate_gcro_noBE[which(replicate_gcro_noBE$Party == "AFRICAN NATIONAL CONGRESS"), ]
#   replicate_gcro_ANC_muns <- replicate_gcro_ANC[-which(replicate_gcro_ANC$Municipality == "GT422 - Midvaal [Meyerton]"), ] # only non-ANC municipality in GT
#   replicate_gcro_DA <- replicate_gcro_noBE[which(replicate_gcro_noBE$Party == "DEMOCRATIC ALLIANCE"), ]
#   replicate_gcro_DA_noMV <- replicate_gcro_DA[-which(replicate_gcro_DA$Municipality == "GT422 - Midvaal [Meyerton]"), ]
#   replicate_gcro_nfs2 <- replicate_gcro_noBE[-which(replicate_gcro_noBE$water_source_13 > 0.9 & replicate_gcro_noBE$elec_13 > 0.9 & replicate_gcro_noBE$flush_toilet_13 > 0.9 & replicate_gcro_noBE$waste_disposal_13 > 0.9), ]
#   replicate_gcro_nfs2_ANC <- replicate_gcro_nfs2[which(replicate_gcro_nfs2$Party.Name == "AFRICAN NATIONAL CONGRESS"), ]
#   
#   boot.model <- glm(formula = renom ~ obj_sp_delta_index_1315 + local_satisfied_15 + 
#                       subj_sp_index15 + formal_housing_delta_1315 + unemp_delta_1315_reverse + 
#                       high_income_avg + post_secondary_avg + black_avg + ce_index15 + 
#                       logpop2011 + years_incumbent + ever_switchedparty 
#                     + 
#                       I(Municipality), family = binomial(link = "logit"), data = replicate_gcro_noBE)
#   
#   results[, b] <- boot.model$coefficients
#   
#   print(b)
# }
# 
# quantile(results[2,], c(0.025, 0.5, 0.975)) #objective service provision change
# quantile(results[3,], c(0.025, 0.5, 0.975)) #councillor satisfaction
# #main results hold up
# 
# # pdf(file="./nominations_paper/bootstrapestimates.pdf")
# pdf("./Outputs/nominations_paper/Replication2/bootstrapestimates.pdf")
# par(mfrow=c(2,1))
# par(mar=c(4,4,2,2))
# plot(density(results[2,]), 
#      main = "Bootstrap estimated coefficients: Service Coverage Change")
# abline(v=quantile(results[2,], c(0.025, 0.5, 0.975))[1], lty=2)
# abline(v=quantile(results[2,], c(0.025, 0.5, 0.975))[3], lty=2)
# abline(v=quantile(results[2,], c(0.025, 0.5, 0.975))[2]) 
# abline(v=0, lwd=2)
# 
# plot(density(results[3,]), 
#      main = "Bootstrap estimated coefficients: Satisfaction with Councillor")
# abline(v=quantile(results[3,], c(0.025, 0.5, 0.975))[1], lty=2)
# abline(v=quantile(results[3,], c(0.025, 0.5, 0.975))[3], lty=2)
# abline(v=quantile(results[3,], c(0.025, 0.5, 0.975))[2]) 
# abline(v=0, lwd=2)
# dev.off()
# 

## Appendix Figure 4: Predicting 2016 win probability using 2011 vote margin for ANC ##

# WARD LEVEL

# Merge 2016 boundary IDs in
load("max_overlap_2011_to_2016.RData")
gcro_noBE <- merge(gcro_noBE, max_overlap_2011_to_2016, by.x = "wardnumber", by.y = "WardID2011", all.x = T)
# Merge 2016 vote margins in
margin2016 <- read.dta13("LGE2016wardshares.dta")
# subset to ANC
margin2016wardANC <- margin2016[margin2016$ballottype == "Ward" & margin2016$party == "AFRICAN NATIONAL CONGRESS", ]
margin2016wardANC$win <- ifelse(margin2016wardANC$voteshareward > 0.5, 1, 0)
margin2016wardANC$wardnumber <- unlist(lapply(strsplit(margin2016wardANC$ward, "Ward "), '[', 2))
# merge
gcro_noBE <- merge(gcro_noBE, margin2016wardANC, by.x = "WardID2016max", by.y = "wardnumber", all.x = T)
# subset to 2011 wards that have more than an 80% overlap with a 2016 ward
gcro_noBE_overlap80 <- gcro_noBE[gcro_noBE$overlap_pct > 0.8, ]
gcro_noBE_overlap80_anc <- gcro_noBE_overlap80[gcro_noBE_overlap80$Party.Name == "AFRICAN NATIONAL CONGRESS", ]

# MUNICIPAL LEVEL

# Merge 2016 vote margins in
margin2016 <- read.dta13("LGE2016wardshares.dta")
# Load municipal-level covariates
muncovars <- read.csv("muncovars.csv")
# Subset to ANC vote share
margin2016_anc_pr <- margin2016[margin2016$party == "AFRICAN NATIONAL CONGRESS" & margin2016$ballottype == "PR", ]
margin2016_anc_pr_sub <- margin2016_anc_pr[, c("municipality", "votesharemuni")]
margin2016_anc_pr_muni <- unique(margin2016_anc_pr_sub)
# If ANC won in 2016
margin2016_anc_pr_muni$anc_win <- ifelse(margin2016_anc_pr_muni$votesharemuni > 0.5, 1, 0)
margin2016_anc_pr_muni$cat_b <- unlist(lapply(strsplit(margin2016_anc_pr_muni$municipality, " - "), '[', 1))
# Merge into municipal-level covariates
muncovars <- merge(muncovars, margin2016_anc_pr_muni, by = "cat_b", all.x = T)
# Subset to municipalities controlled by the ANC in 2011
muncovars_anc <- muncovars[muncovars$winner == "AFRICAN NATIONAL CONGRESS", ]

# Create plots
pdf("vote_prediction_gauteng_paper.pdf")
par(mfrow = c(2, 1))
# reset plot margins to default
par(mar = c(5.1,4.1,4.1,2.1))
# ward
margincuts <- data.frame(start = seq(from = 0, to = 0.90, by = 0.05),
                         end = seq(from = 0.05, to = 0.95, by = 0.05))
winprob <- c()
bootci.low <- c()
bootci.hi <- c()
for(i in 1:nrow(margincuts)){
  mysub <- gcro_noBE_overlap80_anc[gcro_noBE_overlap80_anc$win_margin2011 > margincuts[i, 1] & gcro_noBE_overlap80_anc$win_margin2011 <= margincuts[i, 2], ]
  winprob[i] <- mean(mysub$win, na.rm = T)
  #bootstrap
  # resample with replacement mysub win and take the mean
  rowsampmeans <- c()
  for(j in 1:1000){
    rowsamp <- sample(mysub$win, replace = T)
    rowsampmeans[j] <- mean(rowsamp, na.rm = T)
  }
  bootci.low[i] <- quantile(rowsampmeans, 0.025)
  bootci.hi[i] <- quantile(rowsampmeans, 0.975)
}
par(mfrow = c(1, 1))
plot(x = seq(from = 0.025, to = 0.925, by = 0.05),
     y = winprob,
     xaxt = 'n',
     xlab = "2011 Margin",
     ylab = "Pr(2016 Win)",
     #main = "2011 ANC Vote Margin and 2016 ANC Win Probability \n ANC-Controlled Gauteng Wards",
     cex.main = 0.8,
     col = "gray")
segments(x0 = seq(from = 0.025, to = 0.925, by = 0.05),
         x1 = seq(from = 0.025, to = 0.925, by = 0.05),
         y0 = bootci.low,
         y1 = bootci.hi,
         col = "gray")
axis(side = 1,
     at = seq(from = 0, to = 0.95, by = 0.05))
rug(gcro_noBE_overlap80_anc$win_margin2011)
# divide into thirds
margincuts2 <- data.frame(start = seq(from = 0, to = 0.66, by = 0.33),
                          end = seq(from = 0.33, to = 0.99, by = 0.33))
winprob2 <- c()
bootci.low2 <- c()
bootci.hi2 <- c()
for(i in 1:nrow(margincuts2)){
  mysub <- gcro_noBE_overlap80_anc[gcro_noBE_overlap80_anc$win_margin2011 > margincuts2[i, 1] & gcro_noBE_overlap80_anc$win_margin2011 <= margincuts2[i, 2], ]
  winprob2[i] <- mean(mysub$win, na.rm = T)
  #bootstrap
  # resample with replacement mysub win and take the mean
  rowsampmeans <- c()
  for(j in 1:1000){
    rowsamp <- sample(mysub$win, replace = T)
    rowsampmeans[j] <- mean(rowsamp, na.rm = T)
  }
  bootci.low2[i] <- quantile(rowsampmeans, 0.025)
  bootci.hi2[i] <- quantile(rowsampmeans, 0.975)
}
points(x = c(0.165, 0.495, 0.825),
       y = winprob2,
       pch = 20)
segments(x0 = c(0.165, 0.495, 0.825),
         x1 = c(0.165, 0.495, 0.825),
         y0 = bootci.low2,
         y1 = bootci.hi2,
         lwd = 3)
abline(v = 0.33, lty = 2)
abline(v = 0.66, lty = 2)
dev.off()

# municipal level
pdf("vote_prediction_allmunis_paper.pdf")
margincuts <- data.frame(start = seq(from = 0, to = 0.90, by = 0.05),
                         end = seq(from = 0.05, to = 0.95, by = 0.05))
winprob <- c()
bootci.low <- c()
bootci.hi <- c()
for(i in 1:nrow(margincuts)){
  mysub <- muncovars_anc[muncovars_anc$win_margin > margincuts[i, 1] & muncovars_anc$win_margin <= margincuts[i, 2], ]
  winprob[i] <- mean(mysub$anc_win, na.rm = T)
  #bootstrap
  # resample with replacement mysub win and take the mean
  rowsampmeans <- c()
  for(j in 1:1000){
    rowsamp <- sample(mysub$anc_win, replace = T)
    rowsampmeans[j] <- mean(rowsamp, na.rm = T)
  }
  bootci.low[i] <- quantile(rowsampmeans, 0.025, na.rm = T)
  bootci.hi[i] <- quantile(rowsampmeans, 0.975, na.rm = T)
}

par(mfrow = c(1, 1))
plot(x = seq(from = 0.025, to = 0.925, by = 0.05),
     y = winprob,
     xaxt = 'n',
     ylim = c(0, 1),
     #pch = 20,
     xlab = "2011 Margin",
     ylab = "Pr(2016 Win)",
     #main = "2011 ANC Vote Margin and 2016 ANC Win Probability \n ANC-Controlled Municipalities",
     cex.main = 0.8,
     col = "gray")
segments(x0 = seq(from = 0.025, to = 0.925, by = 0.05),
         x1 = seq(from = 0.025, to = 0.925, by = 0.05),
         y0 = bootci.low,
         y1 = bootci.hi,
         col = "gray")
axis(side = 1,
     at = seq(from = 0, to = 0.95, by = 0.05))
rug(muncovars_anc$win_margin)
# divide into thirds
margincuts2 <- data.frame(start = seq(from = 0, to = 0.66, by = 0.33),
                          end = seq(from = 0.33, to = 0.99, by = 0.33))
winprob2 <- c()
bootci.low2 <- c()
bootci.hi2 <- c()
for(i in 1:nrow(margincuts2)){
  mysub <- muncovars_anc[muncovars_anc$win_margin > margincuts2[i, 1] & muncovars_anc$win_margin <= margincuts2[i, 2], ]
  winprob2[i] <- mean(mysub$anc_win, na.rm = T)
  #bootstrap
  # resample with replacement mysub win and take the mean
  rowsampmeans <- c()
  for(j in 1:1000){
    rowsamp <- sample(mysub$anc_win, replace = T)
    rowsampmeans[j] <- mean(rowsamp, na.rm = T)
  }
  bootci.low2[i] <- quantile(rowsampmeans, 0.025, na.rm = T)
  bootci.hi2[i] <- quantile(rowsampmeans, 0.975, na.rm = T)
}
points(x = c(0.165, 0.495, 0.825),
       y = winprob2,
       pch = 20)
segments(x0 = c(0.165, 0.495, 0.825),
         x1 = c(0.165, 0.495, 0.825),
         y0 = bootci.low2,
         y1 = bootci.hi2,
         lwd = 3)
abline(v = 0.33, lty = 2)
abline(v = 0.66, lty = 2)
dev.off()
